The role of osmorespiratory compromise in hypoxia tolerance of the purportedly oxyconforming teleost Galaxias maculatus

1 📕 READ ME

🚧 UNPUBLISHED AND ONGOING example 🚧 The data in this script are currently unpublished, and the analysis is ongoing. Some sections may be incomplete. This script has been made publicly available to facilitate sharing with collaborators and colleagues who are interested in the progress of the project and to ensure transparency.

SUMMARY This R code estimates the relationship between oxygen consumption (MO₂) and ambient oxygen partial pressure (PO₂) in the Common Galaxias (Galaxias maculatus). It also estimates the critical partial pressure of oxygen for aerobic metabolism (Pcrit), commonly defined as the threshold below which the oxygen consumption rate can no longer be sustained. The associated article is titled “The role of osmorespiratory compromise in hypoxia tolerance of the purportedly oxyconforming teleost Galaxias maculatus.” If you are reading the HTML version of this script, click the Code button in the top right to download the .Rmd file.

AIM The article aims to test whether Galaxias maculatus can maintain oxygen consumption (MO₂) as ambient PO₂ declines and, if so, at what level it reaches the critical partial pressure of oxygen for aerobic metabolism (Pcrit).

AUTHORS Timothy D. Clark ^[a] Luis L. Kuchenmüller ^[a] Elizabeth C. Hoots ^[a] Maryane Gradito ^[a] Jake M. Martin ^[a,b] In no particular order

AFFILIATIONS [a] School of Life and Environmental Sciences, Deakin University, Geelong, VIC, Australia
[b] School of Biological Sciences, Monash University, Clayton, VIC, Australia

CONTRIBUTOR ROLES 🚧 To be added 🚧 Based on the Contributor Roles Taxonomy (CRediT)

DISCLAIMER I (Jake Martin) am dyslexic. I have made an effort to review the script for grammatical errors, but some will likely remain. I apologise. Please feel free to reach out using the contact details below if anything is unclear.

2 📧 Contact

Jake M. Martin

📧 Email: jake.martin@deakin.edu.au

📧 Alt Email: jake.martin.research@gmail.com

🌐 Web: jakemartin.org

🐙 GitHub: JakeMartinResearch

3 📑 Sharing/accessing and citing

1. Licenses/restrictions placed on the data: CC-BY 4.0

2. Link to the associated publication:
   🚧 To be added 🚧

3. Recommended citation for this data:
   🚧 To be added 🚧

4 📦 Required packages and knit settings

A list of the required packages to run the script.

# ---- Install pacman if it's not already installed ----
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman")

# ---- List of required packages ----
pkgs <- c(
  # ----- Data Visualisation -----
  "ggthemes", "bayesplot", "gt", "gtsummary", "plotly", "qqplotr", "gridExtra",
  
  # ----- Tidy Data and Wrangling -----
  "tidyverse", "janitor", "readxl", "broom.mixed", "data.table", "hms", "devtools",
  "mclust",
  
  # ----- Modelling and Statistical Analysis -----
  "brms", "rstan", "marginaleffects", "performance", "emmeans",
  "tidybayes", "respirometry", "future"
)

# ---- Install and load all packages using pacman ----
suppressPackageStartupMessages(
  pacman::p_load(char = pkgs, install = TRUE)
)

#kniter seetting
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE, # no warnings
cache = TRUE,# Cacheing to save time when kniting
tidy = TRUE,
fig.align = "center" 
)

5 🔧 Custom functions

Here are some custom functions used within this script.

bayes_incremental_regression_by_id(): A custom function to build Bayesian incremental regressions. It is designed to run a list of subgroup models (IDs) in parallel using 4 cores. The function uses brm() with a Gaussian error distribution.

Use:
The function accepts the following arguments:
- id_i: A grouping factor or ID used to filter the data for each regression. If none is provided, the function uses the entire dataset.
- id_name: The column name (as a character string) corresponding to the grouping factor in the data frame. If not provided, the function uses all data.
- data: The data frame containing the variables for analysis.
- predictor: The predictor variable of interest.
- response: The response variable of interest.
- seed_number: A random seed value for model reproducibility.
- save_models: A logical argument indicating whether to save the model outputs (TRUE or FALSE).
- mod_output_wd: The output directory where model .rds files should be saved (used only if save_models = TRUE).

# Instead we could also use a distributional regression approach, by
# specifically modelling the variance by DO (e.g. sigma ~ DO). Weighting may
# not be required in this case, I don't think higher density of values in a
# given space will effect Bayesian estimates like it does in frequentist
# models. See discourse https://discourse.mc-stan.org/t/weights-in-brm/4278

bayes_incremental_regression_by_id <- function(id_i, id_name, data, predictor, response,
    seed_number, save_models, mod_output_wd) {
    # Initiate an empty list to store models
    models <- list()

    # Check if id_name is missing, NULL, or blank, and assign NA if so
    if (missing(id_name) || is.null(id_name) || id_name == "") {
        id_name <- NA
    }

    # Check if id_i is missing, NULL, or blank, and assign NA if so
    if (missing(id_i) || is.null(id_i) || id_i == "") {
        id_name <- NA
    }

    # Filter data for the current ID if id_name is given as a factor or
    # character and id_i is defined
    df_i <- data %>%
        dplyr::filter(if (!is.na(id_i) && (is.factor(data[[id_name]]) || is.character(data[[id_name]]))) {
            !!rlang::sym(id_name) == id_i
        } else {
            TRUE
        })

    # Dynamically create formulas
    formula_lm_0 <- reformulate("1", response)
    formula_lm_1 <- reformulate(predictor, response)
    formula_lm_2 <- reformulate(sprintf("poly(%s, 2)", predictor), response)
    formula_lm_3 <- reformulate(sprintf("poly(%s, 3)", predictor), response)

    # Fit and store models in the list
    models[[paste0(id_i, "_lm_0")]] <- brm(bf(formula_lm_0, family = gaussian()),
        data = df_i, cores = 4, seed = seed_number, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_0"))

    models[[paste0(id_i, "_lm_1")]] <- brm(bf(formula_lm_1, family = gaussian()),
        data = df_i, cores = 4, seed = seed_number, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_1"))

    models[[paste0(id_i, "_lm_2")]] <- brm(bf(formula_lm_2, family = gaussian()),
        data = df_i, cores = 4, seed = seed_number, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_2"))

    models[[paste0(id_i, "_lm_3")]] <- brm(bf(formula_lm_3, family = gaussian()),
        data = df_i, cores = 4, seed = 143019, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_3"))

    # Return the list of models for the current ID
    return(models)
}

load_rds9(): A custom function to load all rds models in a directory and store in a list. Takes a directory with .rds files

load_rds <- function(model_dw) {
    # List all .rds files in the directory
    model_file_list <- list.files(path = model_dw, pattern = "\\.rds$", full.names = TRUE)

    # Strip extensions and get names
    model_names <- tools::file_path_sans_ext(basename(model_file_list))

    # Read all .rds files into a named list
    model_store_list <- setNames(lapply(model_file_list, readRDS), model_names)

    return(model_store_list)
}

incremental_regression_bayes_fits(): A custom function for pulling model fits, loo and r2 using loo() and bayes_R2(), respectively. Takes a list of brm models.

# Define Function to Process the data for each ID
incremental_regression_bayes_fits <- function(models) {
  
  loo_results_list <- list()
  
  # Iterate over the names of the models
  for (mod_name in names(models)) {
    # Extract the model
    mod_i <- models[[mod_name]]
    
    # Compute LOO results
    mod_loo_results_i <- loo::loo(mod_i)
    
    # Extract relevant LOO metrics
    elpd_loo_i <- mod_loo_results_i$elpd_loo
    p_loo_i <- mod_loo_results_i$p_loo
    looic_i <- mod_loo_results_i$looic
    
    # Create a data frame with metrics
    df_i <- data.frame(
      elpd_loo = elpd_loo_i,
      p_loo = p_loo_i,
      looic = looic_i,
      model = mod_name
    )
    
    est_i <- tidy(mod_i, effects = "fixed", conf.int = TRUE) %>% 
      dplyr::select(term, estimate, conf.low, conf.high) %>% 
      tidyr::pivot_wider(
        names_from = term,                # Use `term` as column names
        values_from = c(estimate, conf.low, conf.high),  # Values to pivot
        names_sep = "_"                   # Add a separator to column names
      )
    
    df_i <- cbind(df_i, est_i)
    
    # Store the data frame in the list
    loo_results_list[[mod_name]] <- df_i
  }
  
  # Combind 
  loo_results_combined <- bind_rows(loo_results_list)
  
  # Get R2 
  r2_results <- map_dfr(models, ~ as.data.frame(bayes_R2(.x)), .id = "model") %>%
    tibble::remove_rownames()
  
  # Combind R2 and loo results 
  model_fit_df <- dplyr::full_join(loo_results_combined, r2_results, by = "model") %>% 
    dplyr::select(model, everything()) %>% 
    dplyr::rename(r2 = Estimate,
                  r2_error = Est.Error,
                  r2_q2.5 = Q2.5,
                  r2_q97.5 = Q97.5) %>% 
    dplyr::mutate(id = sub("_(lm_\\d+)$", "", model),
                  model_type = sub("^.*_(lm_\\d+)$", "\\1", model))
  
  return(model_fit_df)
}

bayes_mod_predictions(): This function extracts the predicted values using fitted() from a list of models and combines them with the original data file used for the model. These are the posterior mean fitted values (i.e. the expected value of the response variable given the predictor variables and the estimated posterior distributions of the parameters) for each observation in the dataset, along with 95% credible intervals.

bayes_mod_predictions <- function(models, original_data) {

    prediction_list <- list()

    for (mod_name in names(models)) {
        # Extract mod
        mod_i <- models[[mod_name]]

        # Make mode predictions
        model_predictions_i <- as.data.frame(fitted(mod_i, summary = TRUE)) %>%
            dplyr::mutate(model = mod_name, id = sub("_(lm_\\d+)$", "", mod_name),
                model_type = sub("^.*_(lm_\\d+)$", "\\1", mod_name)) %>%
            dplyr::rename(pred_lower = Q2.5, pred_upper = Q97.5, predicted = Estimate,
                pred_error = Est.Error) %>%
            dplyr::select(model, everything())

        id_i <- model_predictions_i$id[1]

        original_data_i <- original_data %>%
            dplyr::filter(id == id_i) %>%
            dplyr::select(-id)

        model_predictions_original_i <- cbind(model_predictions_i, original_data_i)

        prediction_list[[mod_name]] <- model_predictions_original_i
    }
    predictions_df <- bind_rows(prediction_list)
    return(predictions_df)
}

calcSMR(): authored by Denis Chabot used to estimate SMR with several different methods Claireaux and Chabot (2016) ^[1]

^[1] Claireaux, G. and Chabot, D. (2016) Responses by fishes to environmental hypoxia: integration through Fry’s concept of aerobic metabolic scope. Journal of Fish Biology https://doi.org/10.1111/jfb.12833

calcSMR = function(Y, q = c(0.1, 0.15, 0.2, 0.25, 0.3), G = 1:4) {
    u = sort(Y)
    the.Mclust <- Mclust(Y, G = G)
    cl <- the.Mclust$classification
    # sometimes, the class containing SMR is not called 1 the following
    # presumes that when class 1 contains > 10% of cases, it contains SMR,
    # otherwise we take class 2
    cl2 <- as.data.frame(table(cl))
    cl2$cl <- as.numeric(levels(cl2$cl))
    valid <- cl2$Freq >= 0.1 * length(time)
    the.cl <- min(cl2$cl[valid])
    left.distr <- Y[the.Mclust$classification == the.cl]
    mlnd = the.Mclust$parameters$mean[the.cl]
    CVmlnd = sd(left.distr)/mlnd * 100
    quant = quantile(Y, q)
    low10 = mean(u[1:10])
    low10pc = mean(u[6:(5 + round(0.1 * (length(u) - 5)))])
    # remove 5 outliers, keep lowest 10% of the rest, average Herrmann & Enders
    # 2000
    return(list(mlnd = mlnd, quant = quant, low10 = low10, low10pc = low10pc, cl = cl,
        CVmlnd = CVmlnd))
}

calcO2crit(): authored by Denis Chabot used to estimate O2crit (Pcript). Claireaux and Chabot (2016) ^[1]

Note: O2 is assumed to be in percentage of dissolved oxygen (DO)

calcO2crit <- function(Data, SMR, lowestMO2 = NA, gapLimit = 4, max.nb.MO2.for.reg = 20) {
    # AUTHOR: Denis Chabot, Institut Maurice-Lamontagne, DFO, Canada first
    # version written in June 2009 last updated in January 2015
    method = "LS_reg"  # will become 'through_origin' if intercept is > 0
    if (is.na(lowestMO2))
        lowestMO2 = quantile(Data$MO2[Data$DO >= 80], p = 0.05)
    # Step 1: identify points where MO2 is proportional to DO
    geqSMR = Data$MO2 >= lowestMO2
    pivotDO = min(Data$DO[geqSMR])
    lethal = Data$DO < pivotDO
    N_under_SMR = sum(lethal)  # points available for regression?
    final_N_under_SMR = lethal  # some points may be removed at Step 4
    lastMO2reg = Data$MO2[Data$DO == pivotDO]  # last MO2 when regulating
    if (N_under_SMR > 1)
        theMod = lm(MO2 ~ DO, data = Data[lethal, ])
    # Step 2, add one or more point at or above SMR 2A, when there are fewer
    # than 3 valid points to calculate a regression
    if (N_under_SMR < 3) {
        missing = 3 - sum(lethal)
        not.lethal = Data$DO[geqSMR]
        DOlimit = max(sort(not.lethal)[1:missing])  # highest DO acceptable
        # to reach a N of 3
        addedPoints = Data$DO <= DOlimit
        lethal = lethal | addedPoints
        theMod = lm(MO2 ~ DO, data = Data[lethal, ])
    }
    # 2B, add pivotDO to the fit when Step 1 yielded 3 or more values?
    if (N_under_SMR >= 3) {
        lethalB = Data$DO <= pivotDO  # has one more value than 'lethal'
        regA = theMod
        regB = lm(MO2 ~ DO, data = Data[lethalB, ])
        large_slope_drop = (coef(regA)[2]/coef(regB)[2]) > 1.1  # arbitrary
        large_DO_gap = (max(Data$DO[lethalB]) - max(Data$DO[lethal])) > gapLimit
        tooSmallMO2 = lastMO2reg < SMR
        if (!large_slope_drop & !large_DO_gap & !tooSmallMO2)
            {
                lethal = lethalB
                theMod = regB
            }  # otherwise we do not accept the additional point
    }
    # Step 3 if the user wants to limit the number of points in the regression
    if (!is.na(max.nb.MO2.for.reg) & sum(lethal) > max.nb.MO2.for.reg) {
        Ranks = rank(Data$DO)
        lethal = Ranks <= max.nb.MO2.for.reg
        theMod = lm(MO2 ~ DO, data = Data[lethal, ])
        final_N_under_SMR = max.nb.MO2.for.reg
    }
    # Step 4
    predMO2 = as.numeric(predict(theMod, data.frame(DO = Data$DO)))
    Data$delta = (Data$MO2 - predMO2)/predMO2 * 100  # residuals set to zero
    # when below pivotDO
    Data$delta[Data$DO < pivotDO | lethal] = 0
    tol = 0  # any positive residual is unacceptable
    HighValues = Data$delta > tol
    Ranks = rank(-1 * Data$delta)
    HighMO2 = HighValues & Ranks == min(Ranks)  # keep largest residual
    if (sum(HighValues) > 0)
        {
            nblethal = sum(lethal)
            Data$W = NA
            Data$W[lethal] = 1/nblethal
            Data$W[HighMO2] = 1
            theMod = lm(MO2 ~ DO, weight = W, data = Data[lethal | HighMO2, ])
            # This new regression is always an improvement, but there can still
            # be points above the line, so we repeat
            predMO2_2 = as.numeric(predict(theMod, data.frame(DO = Data$DO)))
            Data$delta2 = (Data$MO2 - predMO2_2)/predMO2_2 * 100
            Data$delta2[Data$DO < pivotDO] = 0
            tol = Data$delta2[HighMO2]
            HighValues2 = Data$delta2 > tol
            if (sum(HighValues2) > 0)
                {
                  Ranks2 = rank(-1 * Data$delta2)
                  HighMO2_2 = HighValues2 & Ranks2 == 1  # keep the largest residual
                  nblethal = sum(lethal)
                  Data$W = NA
                  Data$W[lethal] = 1/nblethal
                  Data$W[HighMO2_2] = 1
                  theMod2 = lm(MO2 ~ DO, weight = W, data = Data[lethal | HighMO2_2,
                    ])
                  # is new slope steeper than the old one?
                  if (theMod2$coef[2] > theMod$coef[2]) {
                    theMod = theMod2
                    HighMO2 = HighMO2_2
                  }
                }  # end second search for high value
        }  # end first search for high value
    Coef = coefficients(theMod)
    # Step 5, check for positive intercept
    AboveOrigin = Coef[1] > 0
    # if it is, we use a regression that goes through the origin
    if (AboveOrigin) {
        theMod = lm(MO2 ~ DO - 1, data = Data[lethal, ])
        Coef = c(0, coefficients(theMod))  # need to add the intercept (0)
        # manually to have a pair of coefficients
        method = "through_origin"
        HighMO2 = rep(FALSE, nrow(Data))  # did not use the additional value
        # from Step 4
    }
    po2crit = as.numeric(round((SMR - Coef[1])/Coef[2], 1))
    sum_mod = summary(theMod)
    anov_mod = anova(theMod)
    O2CRIT = list(o2crit = po2crit, SMR = SMR, Nb_MO2_conforming = N_under_SMR, Nb_MO2_conf_used = final_N_under_SMR,
        High_MO2_required = sum(HighMO2) == 1, origData = Data, Method = method,
        mod = theMod, r2 = sum_mod$r.squared, P = anov_mod$"Pr(>F)", lethalPoints = which(lethal),
        AddedPoints = which(HighMO2))
}  # end function

plotO2crit(): authored by Denis Chabot, used to plot the modes used for the calcO2crit() function. Claireaux and Chabot (2016) ^[1]

Note: I added abline(h=lowestMO2, col=“pink”) so that I could visualise the lowestMO2 position

plotO2crit <- function(o2critobj, plotID = "", Xlab = "Dissolved oxygen (% sat.)",
    Ylab = "dotitalumol", smr.cex = 0.9, o2crit.cex = 0.9, plotID.cex = 1.2, Transparency = T,
    ...) {
    # AUTHOR: Denis Chabot, Institut Maurice-Lamontagne, DFO, Canada first
    # version written in June 2009 last updated 2015-02-09 for R plotting
    # devices that do not support transparency (e.g., postscript), set
    # Transparency to FALSE
    smr = o2critobj$SMR
    if (Ylab %in% c("dotitalumol", "italumol", "dotumol", "umol", "dotitalmg", "italmg",
        "dotmg", "mg")) {
        switch(Ylab, dotitalumol = {
            mo2.lab = expression(paste(italic(dot(M))[O[2]], " (", mu, "mol ", O[2],
                " ", min^-1, " ", kg^-1, ")"))
        }, italumol = {
            mo2.lab = expression(paste(italic(M)[O[2]], " (", mu, "mol ", O[2], " ",
                min^-1, " ", kg^-1, ")"))
        }, dotumol = {
            mo2.lab = expression(paste(dot(M)[O[2]], " (", mu, "mol ", O[2], " ",
                min^-1, " ", kg^-1, ")"))
        }, umol = {
            mo2.lab = expression(paste(M[O[2]], " (", mu, "mol ", O[2], " ", min^-1,
                " ", kg^-1, ")"))
        }, dotitalmg = {
            mo2.lab = expression(paste(italic(dot(M))[O[2]], " (mg ", O[2], " ",
                h^-1, " ", kg^-1, ")"))
        }, italmg = {
            mo2.lab = expression(paste(italic(M)[O[2]], " (mg ", O[2], " ", h^-1,
                " ", kg^-1, ")"))
        }, dotmg = {
            mo2.lab = expression(paste(dot(M)[O[2]], " (mg ", O[2], " ", h^-1, " ",
                kg^-1, ")"))
        }, mg = {
            mo2.lab = expression(paste(M[O[2]], " (mg ", O[2], " ", h^-1, " ", kg^-1,
                ")"))
        })
    } else mo2.lab = Ylab
    if (Transparency) {
        Col = c(rgb(0, 0, 0, 0.7), "red", "orange")
    } else {
        Col = c(grey(0.3), "red", "orange")
    }
    Data = o2critobj$origData
    lowestMO2 = quantile(Data$MO2[Data$DO >= 80], p = 0.05)  # I added this
    Data$Color = Col[1]
    Data$Color[o2critobj$lethalPoints] = Col[2]
    Data$Color[o2critobj$AddedPoints] = Col[3]
    # ordinary LS regression without added points: blue line, red symbols
    # ordinary LS regression with added points: blue line, red & orange symbols
    # regression through origin: green dotted line, red symbols
    line.color = ifelse(o2critobj$Method == "LS_reg", "blue", "darkgreen")
    line.type = ifelse(o2critobj$Method == "LS_reg", 1, 3)
    limX = c(0, max(Data$DO))
    limY = c(0, max(Data$MO2))
    plot(MO2 ~ DO, data = Data, xlim = limX, ylim = limY, col = Data$Color, xlab = Xlab,
        ylab = mo2.lab, ...)
    coord <- par("usr")
    if (plotID != "") {
        text(0, coord[4], plotID, cex = plotID.cex, adj = c(0, 1.2))
    }
    abline(h = lowestMO2, col = "pink")  # I added this
    abline(h = smr, col = "orange")
    text(coord[1], smr, "SMR", adj = c(-0.1, 1.3), cex = smr.cex)
    text(coord[1], smr, round(smr, 1), adj = c(-0.1, -0.3), cex = smr.cex)
    if (!is.na(o2critobj$o2crit)) {
        abline(o2critobj$mod, col = line.color, lty = line.type)
        segments(o2critobj$o2crit, smr, o2critobj$o2crit, coord[3], col = line.color,
            lwd = 1)
        text(x = o2critobj$o2crit, y = 0, o2critobj$o2crit, col = line.color, cex = o2crit.cex,
            adj = c(-0.1, 0.5))
    }
}  # end of function

6 📂 Directories

6.1 Input

📥 input_data_wd: Directory for the metadata

wd <- getwd()
input_data_wd <- paste0(wd, "./input-data")  # creates a variable with the name of the wd we want to use

📥 mod_data_wd: Directory for model output data estimated slopes

mod_data_wd <- paste0(wd, "./mod-data")

6.2 Output

📤 output_fig_wd: this is where we will put the figures

output_fig_wd <- paste0(wd, "./output-fig")
if (!dir.exists("output-fig")) {
    dir.create("output-fig")
}

📤 output_mods_wd: this is where we will put the models

output_mods_wd <- paste0(wd, "./output-mod")
if (!dir.exists("output-mod")) {
    dir.create("output-mod")
}

7 💿 Data

7.1 Slopes (ṀO)

💿 slope_df: We have imported the slopes extracted in LabChart during each phase of the experiment

 setwd(input_data_wd)

data_file <- "labchart-slope-data.xlsx"

# 
# # Get the names of all sheets in the Excel file
sheet_names <- readxl::excel_sheets(data_file)
all_trials_select <- c("start_date", "order", "phase", "cycle", "date", "time")
slope_list <- list()

for (sheet in sheet_names) {

  df <- read_excel(data_file, sheet = sheet) %>% 
  dplyr::rename_with(tolower)
  
a_name <- paste0("a_", tolower(sheet))
a_df <- df %>%
  dplyr::select(starts_with('a'), all_trials_select) %>% 
  dplyr::rename(temp = a_temp) %>% 
  dplyr::mutate(across(starts_with('a'), as.numeric)) %>% 
  pivot_longer(
    cols = starts_with('a'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>%
  dplyr::mutate(respirometer_group = "a") # Add a new column with a fixed value

slope_list[[a_name]]<- a_df

b_name <- paste0("b_", tolower(sheet))
b_df <- df %>% 
  dplyr::select(starts_with('b'), all_trials_select) %>% 
  dplyr::rename(temp = b_temp) %>% 
  dplyr::mutate(across(starts_with('b'), as.numeric)) %>% 
  pivot_longer(
    cols = starts_with('b'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>% 
    dplyr::mutate(respirometer_group = "b")

slope_list[[b_name]] <- b_df

c_name <- paste0("c_", tolower(sheet))
c_df <- df %>% 
  dplyr::select(starts_with('c'), all_trials_select) %>% 
  dplyr::rename(temp = c_temp,
                i_cycle = cycle) %>% 
  dplyr::mutate(across(starts_with('c'), as.numeric)) %>%
  pivot_longer(
    cols = starts_with('c'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>% 
    dplyr::mutate(respirometer_group = "c") %>% 
  dplyr::rename(cycle = i_cycle)

slope_list[[c_name]] <- c_df

d_name <- paste0("d_", tolower(sheet))
d_df <- df %>% 
  dplyr::select(starts_with('d'), all_trials_select) %>% 
  dplyr::rename(temp = d_temp,
                i_date = date) %>% 
  dplyr::mutate(across(starts_with('d'), as.numeric)) %>%
  pivot_longer(
    cols = starts_with('d'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>% 
    dplyr::mutate(respirometer_group = "d") %>% 
  dplyr::rename(date = i_date)

slope_list[[d_name]] <- d_df
}


slope_df <- bind_rows(slope_list) %>% 
  dplyr::mutate(resp_cat_date = paste0(respirometer_group, "_", start_date),
                chamber_n = str_extract(chamber_id, "\\d+"),
                id_prox = paste0(resp_cat_date, "_", chamber_n),
                time_hms = as_hms(time*3600),
                date_chr = format(date, "%d/%m/%Y")
                )

7.2 Fish metadata

💿 metadata: This is the meta data for each chamber

Note: We are also adding volume based on chamber type.

setwd(input_data_wd)

metadata <- read_excel("fish-meta.xlsx", na = "NA") %>%
    dplyr::mutate(id_split = id) %>%
    tidyr::separate(id_split, into = c("respirometer_group", "salinity_group", "start_date",
        "chamber"), sep = "_") %>%
    dplyr::mutate(volume = dplyr::case_when(chamber_type == "L" ~ 0.3, chamber_type ==
        "M_M" ~ 0.105, chamber_type == "M_NM" ~ 0.11, chamber_type == "S" ~ 0.058,
        chamber_type == "SM" ~ 0.075, chamber_type == "D3" ~ 0.055, TRUE ~ NA), id_prox = paste0(respirometer_group,
        "_", start_date, "_", chamber), rel_size = mass/volume)

7.3 Urbina, Glover, and Forster (2012)

💿 urbina_et_al_2012: This is the mean level data extracted from Urbina, Glover, and Forster (2012)^[2] Figure 1a. We used the metaDigitise package in R to extract the data ^[3].

setwd(input_data_wd)

urbina_et_al_2012 <- read_excel("urbina-et-al-2012-fig1a.xlsx")

7.4 O₂crit viusal inspection

💿 o2crit_check: This data frame includes the visual based assessment for O₂crit for all fish included in the analysis (n = 58). The inspections were made using figures genbrated by the plotO2crit() function, and can be viewed in combined_chabot_plots.pdf. The visual assessment was done independently by all authors, and the presence of an O₂crits allocated to ‘yes’, ‘no’, ‘maybe’. If yes, or maybe, an estimate of the dissolved oxygen percentage was given.

Directory to the figures

📁 gmac-mo2-24/ └── 📂 output-fig/
└── 📂 model_chabot/ └── 📄 combined_chabot_plots.pdf

setwd(input_data_wd)

o2crit_check <- read_excel("o2crit-visual-assessment.xlsx", na = "NA")

7.4.1 Combinding metadata

Adding the meta data to the slopes data frame

slope_df_2 <- slope_df %>%
    dplyr::select(-start_date, -respirometer_group) %>%
    left_join(metadata, by = "id_prox") %>%
    dplyr::mutate(light_dark = if_else(time_hms >= as.hms("07:00:00") & time_hms <
        as.hms("19:00:00"), "light", "dark")) %>%
    dplyr::arrange(id)

8 🧹 Data tidy

8.1 Numbers

We have 64 fish with MO₂ data

n <- slope_df_2 %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::distinct(id) %>%
    nrow(.)

paste0("n = ", n)

## [1] "n = 64"

With 48 from the 0 ppt and 48 from 9 ppt groups

slope_df_2 %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::group_by(salinity_group) %>%
    dplyr::reframe(`n total` = length(unique(id))) %>%
    gt() %>%
    cols_label(salinity_group = "Salinity group") %>%
    cols_align(align = "center", columns = everything())

Salinity group	n total
0	33
9	31

8.2 Fish size

Here we calculate the mean length and size of fish used in the experiment.

mass_length <- slope_df_2 %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::group_by(id) %>%
    dplyr::sample_n(1) %>%
    dplyr::ungroup() %>%
    dplyr::reframe(x_mass = round(mean(mass, na.rm = TRUE), 3), min_mass = round(min(mass,
        na.rm = TRUE), 3), max_mass = round(max(mass, na.rm = TRUE), 3), x_length = round(mean(length,
        na.rm = TRUE), 2), min_length = round(min(length, na.rm = TRUE), 2), max_length = round(max(length,
        na.rm = TRUE), 2))

mass_mean <- mass_length %>%
    pull(x_mass)

mass_min <- mass_length %>%
    pull(min_mass)

mass_max <- mass_length %>%
    pull(max_mass)

length_mean <- mass_length %>%
    pull(x_length)

length_min <- mass_length %>%
    pull(min_length)

length_max <- mass_length %>%
    pull(max_length)

paste0("The mean mass of fish was ", mass_mean, " g (range: ", mass_min, "-", mass_max,
    ")", ", and the mean length was ", length_mean, " mm (range: ", length_min, "-",
    length_max, ")")

## [1] "The mean mass of fish was 0.532 g (range: 0.21-1.6), and the mean length was 50.41 mm (range: 40-70)"

8.3 Filtering trials

We will remove 6 trials which had errors. These are as follows:

• a_0_25nov_3 fish died during trial
• b_0_26nov_4 flat lined early
• c_0_22nov_2 chamber was opened early
• c_9_26nov_2 chamber was opened early
• c_9_26nov_4 chamber was opened early
• d_9_27nov_3 sensor was jumpy and end points were hard to confidently ID visually

remove_trial_error <- c("a_0_25nov_3", "b_0_26nov_4", "c_0_22nov_2", "c_9_26nov_2",
    "c_9_26nov_4", "d_9_27nov_3")

slope_df_filter <- slope_df_2 %>%
    dplyr::filter(!(id %in% remove_trial_error))

We now have 58 fish with MO2 data

n <- slope_df_filter %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::distinct(id) %>%
    nrow(.)

paste0("n = ", n)

## [1] "n = 58"

With 30 in the 0 ppt group and 28 in the 9 ppt group

slope_df_filter %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::group_by(salinity_group) %>%
    dplyr::reframe(`n total` = length(unique(id))) %>%
    gt() %>%
    cols_label(salinity_group = "Salinity group") %>%
    cols_align(align = "center", columns = everything())

Salinity group	n total
0	30
9	28

Chamber size breakdown

slope_df_filter %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::mutate(volume = if_else(volume == 0.055, 0.058, volume)) %>%
    dplyr::group_by(volume) %>%
    dplyr::reframe(`n total` = length(unique(id)), median_mass = median(mass, na.rm = TRUE),
        sd_mass = sd(mass, na.rm = TRUE)) %>%
    gt() %>%
    cols_label(volume = "Chamber type", median_mass = "Median mass", sd_mass = "Standard deviation mass") %>%
    cols_align(align = "center", columns = everything())

Chamber type	n total	Median mass	Standard deviation mass
0.058	36	0.43	0.1496554
0.105	21	0.61	0.2561871
0.300	1	1.60	0.0000000

Numbers per trial approach

phases <- c("50c", "100c")

slope_df_filter %>%
    dplyr::filter(chamber_condition == "fish" & phase %in% phases) %>%
    dplyr::group_by(phase) %>%
    dplyr::reframe(`n total` = length(unique(id))) %>%
    gt() %>%
    cols_label(phase = "Type of trial", ) %>%
    cols_align(align = "center", columns = everything())

Type of trial	n total
100c	19
50c	39

8.4 Filtering MO₂ estimates

Here we apply the following filters to the MO₂ data:

• Remove the first 5 SMR cycles (burn in)
• Remove all positive raw slopes
• Remove all MO₂ calculated using less then 60 data points (5 min at 0.2 Hz)
• Remove all MO₂ calculated if O₂ increases in a closed phase (i.e. trial has ended)

Check positive values for MO₂ before removing.

slope_tidy_remove_flush <- slope_df_filter %>%
  dplyr::filter(phase != "smr", chamber_condition == "fish") %>% 
  dplyr::group_by(id) %>%
  dplyr::arrange(id, order) %>%  # Ensure the data is sorted within each group
  dplyr::mutate(o2_diff = if_else(row_number() == 1, 0, o2 - lag(o2)), # Calculate the difference in 'o2'
                o2_diff_cumsum = cumsum(o2_diff > 1)) %>%  # Checks first occurrence and sums 
  dplyr::filter(o2_diff_cumsum == 0) %>%  # Keep rows until the first jump > 1
  dplyr::ungroup() %>%
  dplyr::select(-o2_diff, -o2_diff_cumsum)

postive_slopes <- slope_tidy_remove_flush %>%
  dplyr::filter(chamber_condition == "fish" & mo2corr > 0)

n_postive_slopes <- nrow(postive_slopes)

n_slopes <-  nrow(slope_tidy_remove_flush)

list_postive_all <- postive_slopes %>% 
  dplyr::distinct(id) %>% 
  dplyr::pull(id)

print(paste0("There are ", length(list_postive_all), " fish with postive slopes. These fish are: ", paste(list_postive_all, collapse = ", "), ". For all estimated slopes (n = ", n_slopes, ") ", round((n_postive_slopes/n_slopes)*100,2) , "% were postive (n = ", n_postive_slopes, ")"))

## [1] "There are 13 fish with postive slopes. These fish are: a_0_24nov_1, b_0_24nov_3, b_9_21nov_2, b_9_21nov_3, b_9_22nov_1, b_9_22nov_2, b_9_22nov_3, c_0_21nov_1, c_0_21nov_2, c_0_22nov_4, c_9_24nov_3, d_0_22nov_2, d_9_25nov_3. For all estimated slopes (n = 1105) 1.27% were postive (n = 14)"

Filtering the MO₂ data

cycle_burn <- 0:4

slope_df_filter_1 <- slope_df_filter %>%
  dplyr::filter(!(cycle %in% cycle_burn) & 
                  mo2corr < 0 & 
                  n > 60 &
                  chamber_condition == "fish"
                )
  
# Filter out the end flush
slope_tidy_closed <- slope_df_filter_1 %>%
  dplyr::filter(phase != "smr") %>% 
  dplyr::group_by(id) %>%
  dplyr::arrange(id, order) %>%  # Ensure the data is sorted within each group
  dplyr::mutate(o2_diff = if_else(row_number() == 1, 0, o2 - lag(o2)), # Calculate the difference in 'o2'
                o2_diff_cumsum = cumsum(o2_diff > 1)) %>%  # Checks first occurrence and sums 
  dplyr::filter(o2_diff_cumsum == 0) %>%  # Keep rows until the first jump > 1
  dplyr::ungroup() %>%
  dplyr::select(-o2_diff, -o2_diff_cumsum)

slope_tidy_smr <- slope_df_filter_1 %>% 
  dplyr::filter(phase == "smr")

slope_df_filter_2 <- rbind(slope_tidy_smr, slope_tidy_closed) %>% 
  dplyr::arrange(id, order)

8.5 Calculating SMR

We will calculate SMR using calcSMR function by Chabot, Steffensen and Farrell (2016)^[1]. Specifically, we use mean of the lowest normal distribution (MLND) where CVmlnd < 5.4 and the mean of the lower 20% quantile (q0.2) were CVmlnd > 5.4. If CVmlnd is not calculated we have used q0.2.

labchart_chabot_smr <- slope_df_filter_2 %>%
    dplyr::filter(phase == "smr")

# Extract distinct IDs
ids <- labchart_chabot_smr %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

# Initialise an empty list to store SMR data
smr_list <- list()

# Process each ID
for (id_i in ids) {
    tryCatch({
        # Filter the data for the current ID
        df_i <- labchart_chabot_smr %>%
            dplyr::filter(id == id_i) %>%
            dplyr::mutate(abs_mo2corr = abs(mo2corr))

        # Calculate SMR results
        calcSMR_results <- calcSMR(df_i$abs_mo2corr)
        CVmlnd_i <- calcSMR_results$CVmlnd
        quant_i <- calcSMR_results$quant %>%
            as_tibble()
        quant_20per_i <- quant_i$value[3]
        mlnd_i <- calcSMR_results$mlnd
        smr_value <- if_else(CVmlnd_i < 5.4, mlnd_i, quant_20per_i)
        smr_type <- if_else(CVmlnd_i < 5.4, "mlnd", "quant_20per")
        smr_value <- if_else(is.na(smr_value), quant_20per_i, smr_value)
        smr_type <- if_else(is.na(smr_type), "quant_20per", smr_type)

        # Create a data frame for the current ID
        smr_df_i <- tibble::tibble(id = id_i, smr = smr_value, smr_est = smr_type)

    }, error = function(e) {
        # Handle errors by assigning NA values
        smr_df_i <- tibble::tibble(id = id_i, smr = NA, smr_est = NA)
    })

    # Append to the list
    smr_list[[id_i]] <- smr_df_i
}

# Combine all individual SMR data frames into one
smr_df <- bind_rows(smr_list) %>%
    dplyr::rename(smr_chabot = smr, smr_chabot_method = smr_est)

slope_df_filter_3 <- slope_df_filter_2 %>%
    dplyr::left_join(., smr_df, by = "id")

8.6 Transforming MO₂

Here we are transforming the MO₂ units. The resulting values are as follows:

• MO2 is absolute value of the background and leak corrected MO₂ slope from Labchart (mo2corr) times the net volume of the chamber (volume - fish mass), × 60, × 60, to achieve MO₂ as mg^-1 O₂ h^-1

• MO2_g is MO2 divided by fish mass to achieve MO₂ as mg^-1 O₂ g^-1 h^-1 (i.e. mass standardised)

• SMR absolute value of the SMR estimates using methods described by Chabot, Steffensen and Farrell (2016)^[1] times the net volume of the chamber (volume - fish mass), × 60, × 60, to achieve SMR as mg^-1 O₂ h^-1)

• SMR_g is SMR divided by fish mass to achieve SMR as mg^-1 O₂ g^-1 h^-1 (i.e. mass standardised)

• DO is dissolved oxygen percentage calculated from O₂ values (mg^-1 L^-1) using the recorded temperature, salinity, and a constant atmospheric pressure (Pa; 1013.25)

• o2_kpa is the O₂ concentration in kilopascal (kpa). This is used to make a comparative figure only.

# Combine back into one data frame
slope_tidy <- slope_df_filter_3 %>% 
    dplyr::mutate(DO = conv_o2(
                   o2 = o2,
                   from = "mg_per_l",
                   to = "percent_a.s.",
                   temp = temp, #C
                   sal = measured_salinity,
                   atm_pres = 1013.25),
                  o2_kpa = conv_o2(
                   o2 = o2,
                   from = "mg_per_l",
                   to = "kPa",
                   temp = temp, #C
                   sal = measured_salinity,
                   atm_pres = 1013.25),
                  net_volume = volume - mass/1000,
                  MO2 = abs(mo2corr)*net_volume*60*60,
                  MO2_g = MO2/mass,
                  SMR = abs(smr_chabot)*net_volume*60*60,
                  SMR_g = SMR/mass
                  )

Trial lengths

slope_tidy %>%
    dplyr::filter(phase != "smr") %>%
    dplyr::mutate(datetime = as.POSIXct(paste(date, time_hms), format = "%Y-%m-%d %H:%M:%S")) %>%
    dplyr::group_by(id) %>%
    dplyr::reframe(phase = phase[1], start_time = min(datetime, na.rm = TRUE), end_time = max(datetime,
        na.rm = TRUE), time_diff = difftime(end_time, start_time, units = "hours")  # Change units as needed
) %>%
    dplyr::group_by(phase) %>%
    dplyr::reframe(`Minimum duration` = round(min(time_diff), 2), `Maximum duration` = round(max(time_diff),
        2), `Median duration` = round(median(time_diff), 2)) %>%
    dplyr::mutate(phase = if_else(phase == "100c", "Closed at 100% air saturation",
        "Closed at 75% and 50% air saturation")) %>%
    dplyr::rename(Approach = phase) %>%
    gt()

Approach	Minimum duration	Maximum duration	Median duration
Closed at 100% air saturation	3	13	7
Closed at 75% and 50% air saturation	2.91	8.02	5.96

9 📊 Exploratory Data Analysis (EDA)

9.1 O₂ vs MO₂

9.1.1 Figure S1

This interactive was used to identify any outliers, or potential errors.

lm_lines <- slope_tidy %>%
  group_by(id) %>%
  summarise(
    DO_seq = list(seq(min(DO), max(DO), length.out = 100)),  # Generate a sequence of DO values
    MO2_pred = list(predict(lm(MO2_g ~ DO, data = cur_data()), newdata = data.frame(DO = seq(min(DO), max(DO), length.out = 100))))
  ) %>%
  unnest(c(DO_seq, MO2_pred))  # Expand lists into rows

# Create scatter plot with markers for each fish
p <- plot_ly(
  data = slope_tidy,
  x = ~DO,
  y = ~MO2_g,
  type = "scatter",
  mode = "markers",
  color = ~id,  # Colour points by fish ID
  marker = list(opacity = 0.6),
  name = ~id
)

# Add regression lines for each fish
p <- p %>%
  add_trace(
    data = lm_lines,
    x = ~DO_seq,
    y = ~MO2_pred,
    type = "scatter",
    mode = "lines",
    color = ~id,  # Ensure each line matches its corresponding fish
    line = list(width = 1, dash = "solid"),
    showlegend = FALSE  # Avoid cluttering the legend
  )

# Final layout
p <- p %>%
  layout(
    title = "MO<sub>2</sub> vs Dissolved Oxygen with individual linear regressions",
    xaxis = list(title = "Dissolved Oxygen (%)"),
    yaxis = list(title = "MO<sub>2</sub> (mg<sup>-1</sup> O<sub>2</sub> g<sup>-1</sup> h<sup>-1</sup>)"),
    showlegend = FALSE
  )

# Display plot
p

Figure S1: interactive plot of metabolic rate measurements (MO₂; mg O₂ g^-1h^-1) by dissolved oxygen percentage (DO) for all fish, including all estimates during the SMR phase (i.e. intermittent phase). Individual linear regression were fitted for visual reference, and do not represent the best fitting regression.

9.1.2 Figure S2

Looking at the difference responses in the two salinity groups.

salinity_summary <- slope_tidy %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(n = length(unique(id)))

slope_tidy %>% 
  ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
  geom_point(show.legend = FALSE) +
  geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
  geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = TRUE, colour = "red", linetype = "dashed") +
  theme_clean() +
  facet_wrap(~salinity_group) +
  labs(
    subtitle = "mo2 vs o2 by salinity treatment",
    x = "Dissolved oxygen percentage (DO)",
    y = "MO2 (O2 mg/g/h)"
  ) +
  geom_text(data = salinity_summary,
            aes(x = -Inf, y = Inf, label = paste0("italic(n) == ", n)),
            hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, parse = TRUE)

Figure S2: Metabolic rate measurements (MO₂; mg O₂ g^-1h^-1) by dissolved oxygen percentage (DO) for fish from the two salinity treatments, including all estimates during the SMR phase (i.e. intermittent phase). Individual linear regression were fitted for visual reference, and do not represent the best fitting regression. The solid red and dashed red line represents the global trend in each of the treatments, both a linear (solid red) and non-linear (dashed red; Generalized Additive Model fitted with geom_smooth).

9.1.3 Figure S3

A plot to look at the different chamber types

chamber_summary <- slope_tidy %>% 
  dplyr::group_by(chamber_type) %>% 
  dplyr::reframe(n = length(unique(id)))

slope_tidy %>% 
  ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
  geom_point(show.legend = FALSE) +
  geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
  geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
  geom_smooth(se = TRUE, colour = "red", linetype = "dashed") +
  theme_clean() +
  facet_wrap(~chamber_type, scale = "free") +
  labs(
    subtitle = "mo2 vs o2 by chamber type",
    x = "Dissolved oxygen percentage (DO)",
    y = "MO2 (mg O2 g/h)"
  ) +
  geom_text(data = chamber_summary,
            aes(x = -Inf, y = Inf, label = paste0("italic(n) == ", n)),
            hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, parse = TRUE)

Figure S3: Metabolic rate measurements (MO₂; mg O₂ g^-1h^-1) by dissolved oxygen percentage (DO) for fish tested in the 4 different chamber types, including all estimates during the SMR phase (i.e. intermittent phase). Individual linear regression were fitted for visual reference, and do not represent the best fitting regression. The solid red and dashed red line represents the global trend in each of the treatments, both a linear (solid red) and non-linear (dashed red; Generalized Additive Model fitted with geom_smooth).

9.2 Routine MO₂

Making an SMR phase only data frame

slope_tidy_smr <- slope_tidy %>%
    dplyr::filter(phase == "smr")

9.2.1 Figure S4

Routine MO₂ by salinity

mean_mo2_salinity <- slope_tidy_smr %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(mean_mo2 = mean(MO2, na.rm = TRUE))

fig_i <- ggplot() +
    geom_violin(data = slope_tidy_smr, aes(x = salinity_group, y = MO2, fill = salinity_group), color = NA, alpha = 0.3) +
  geom_jitter(data = slope_tidy_smr, aes(x = salinity_group, y = MO2, fill = salinity_group),
                       shape = 21, size = 2, color = "black", alpha = 0.2) +
  geom_boxplot(data = slope_tidy_smr, aes(x = salinity_group, y = MO2, fill = salinity_group),
                        size = 1, alpha = 0.5, outlier.shape = NA, width = 0.3) +
  geom_point(data = mean_mo2_salinity, 
                aes(x = salinity_group, y = mean_mo2, fill = salinity_group), 
                size = 3, alpha = 0.8, colour = "black", stroke = 2) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme_clean() +
  theme(legend.position = "none") +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "Routine MO2 (mg O2 g/h)"
  )

fig_i

Figure S4: Plot of all MO₂ measures during SMR phase by salinity treatment. The small points are the raw observed values, the shaded area behind the points is the a kernel density of the observed data, the box plot shows the median and interquartile range (IQR), and the large point shows the mean.

9.3 SMR

9.3.1 Figure S5

Here’s the same plot but for only the SMR, as estimated with calcSMR function by Chabot, Steffensen and Farrell (2016)^[1]. Specifically, we use mean of the lowest normal distribution (MLND) where CVmlnd < 5.4 and the mean of the lower 20% quantile (q0.2) were CVmlnd > 5.4. If CVmlnd is not calculated we have used q0.2.

smr_only <- slope_tidy_smr %>% 
  dplyr::group_by(id) %>% 
  dplyr::slice(1)

mean_smr_salinity <- smr_only %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(mean_smr = mean(SMR, na.rm = TRUE))

fig_i <- ggplot() +
    geom_violin(data = smr_only, aes(x = salinity_group, y = SMR, fill = salinity_group), color = NA, alpha = 0.3) +
  geom_jitter(data = smr_only, aes(x = salinity_group, y = SMR, fill = salinity_group),
                       shape = 21, size = 2, color = "black", alpha = 0.2) +
  geom_boxplot(data = smr_only, aes(x = salinity_group, y = SMR, fill = salinity_group),
                        size = 1, alpha = 0.5, outlier.shape = NA, width = 0.3) +
  geom_point(data = mean_smr_salinity, 
                aes(x = salinity_group, y = mean_smr, fill = salinity_group), 
                size = 3, alpha = 0.8, colour = "black", stroke = 2) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme_clean() +
  theme(legend.position = "none") +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "MO2 (mg O2 g/h)"
  )

fig_i

Figure S5: The standard metabolic rate (SMR) by salinity treatment. The small points are the raw observed values, the shaded area behind the points is the a kernel density of the observed data, the box plot shows the median and interquartile range (IQR), and the large point shows the mean.

9.4 Individual O₂, MO₂, and SMR

Here we will plot the individual relationship between O₂, MO₂, and SMR for all fish

Create output directory if needed

Loop through each fish ID to create a plot, save these to a single PDF file, this is called combined_slopes.pdf

# # Define consistent colors for all phases phase_colors <- c( '50c' =
# '#e56b6f', '75c' = '#b56576', '100c' = '#6d597a', 'smr' = '#355070' ) ids <-
# slope_tidy %>% dplyr::distinct(id) %>% pull(id) %>% as.list() MO2_plot_list
# <- list() # 1) Open the PDF device once pdf( file =
# file.path(output_fig_slopes_wd, 'combined_slopes.pdf'), width = 8, height = 6
# ) # 2) Loop over IDs and create each plot for (id_i in ids) { smr <-
# slope_tidy %>% dplyr::filter(id == id_i) %>% dplyr::slice(1) %>%
# dplyr::pull(SMR) plot <- slope_tidy %>% dplyr::filter(id == id_i) %>%
# ggplot(aes(x = o2, y = MO2)) + geom_hline(yintercept = smr, linetype =
# 'dashed', color = 'darkred') + geom_point(aes(colour = phase)) +
# scale_color_manual(values = phase_colors, drop = FALSE) + # Ensures
# consistent colours theme_clean() + labs( subtitle = paste0(id_i, ' slopes'),
# x = 'Mean o2 (mg_per_l)', y = 'abs(mo2) (mg_per_l)' ) # Instead of saving
# each plot separately, just print it print(plot) MO2_plot_list[[id_i]] <- plot
# } # 3) Close the PDF device *after* the loop dev.off()

Figure S7: Metabolic rate measurements (MO₂; mg O₂ g^-1h^-1 by O₂) at each of the four experimental phase (the over night SMR intermittent-flow respirometry phase, closed phases at 50%, 75% or 100% O₂). The estimated SMR value is represented as a dashed red line.

10 🧮 Analysis

10.1 Routine MO2

10.1.1 Scaling predictors

Here we scale our predictors for the model. We have first applied an log transformation to mass, so its directly interpretable as the mass-scaling exponent, just like in metabolic theory.

scale_list <- c("temp", "order", "log_mass")

slope_tidy_smr <- slope_tidy_smr %>%
    dplyr::mutate(log_mass = log(mass), across(all_of(scale_list), ~scale(.x, center = TRUE,
        scale = FALSE), .names = "{.col}_c"), light_dark_c = if_else(light_dark ==
        "light", 0.5, -0.5))

Summary of predictors used in the model below

slope_tidy_smr %>%
    dplyr::mutate(cycle = as.numeric(cycle)) %>%
    dplyr::reframe(temp_range = paste0(round(min(temp, na.rm = TRUE), 2), "–",
        round(max(temp, na.rm = TRUE), 2)), mass_range = paste0(round(min(mass, na.rm = TRUE),
        2), "–", round(max(mass, na.rm = TRUE), 2)), cycle_range = paste0(round(min(cycle,
        na.rm = TRUE), 2), "–", round(max(cycle, na.rm = TRUE), 2)), number_fish = length(unique(id))) %>%
    gt()

temp_range	mass_range	cycle_range	number_fish
13.84–14.38	0.21–1.6	5–27	58

10.1.2 Model structure

Here we will use a Bayesian Generalised Linear Mixed Model (GLMM) with a Gamma distribution and a log link, where the shape parameter (α) is also modelled as a function of predictors. This models MO₂ by salinity during the SMR phase to see if the fish held at different salinities have different RMRs. We have also added a few scaled predictors, that may help describe variation in the data, such as mass (g; 0.21–1.6) temperature (°C; 13.84–14.38), cycle order (5–27), and light/dark cycle (light or dark; light between 07:00:00 and 19:00:00), we also include a random effect for fish id to account for multiple MO₂ measures on each fish (1–58). We allowed the the shape parameter (α) to vary as a function of some of the predictors (e.g. salinity_group, order_z) to improve fit.

mo2_gamma_bf <- bf(MO2 ~ temp_c + order_c + light_dark_c + log_mass_c + salinity_group +
    (1 | id), shape ~ salinity_group + order_c, family = Gamma(link = "log"))

10.1.3 Prior selection

These are the default priors. We will use these.

default_prior <- get_prior(mo2_gamma_bf, data = slope_tidy_smr, family = Gamma(link = "log"))
default_prior %>%
    gt()

prior	class	coef	group	resp	dpar	nlpar	lb	ub	source
	b								default
	b	light_dark_c							default
	b	log_mass_c							default
	b	order_c							default
	b	salinity_group9							default
	b	temp_c							default
student_t(3, -2.7, 2.5)	Intercept								default
student_t(3, 0, 2.5)	sd						0		default
	sd		id						default
	sd	Intercept	id						default
	b				shape				default
	b	order_c			shape				default
	b	salinity_group9			shape				default
student_t(3, 0, 2.5)	Intercept				shape				default

10.1.4 Run model

Here we run the model, I have hashed this out because I have saved the model for quick reloading.

⏭️ Skip this step if you have downloaded the saved model ‘mo2_mod_gamma.rds’

Here we reload the model

mo2_mod_gamma <- readRDS(file = paste0(output_mods_wd, "./mo2_mod_gamma.rds"))

10.1.5 Model diagnostics

Checking model convergence

color_scheme_set("red")
plot(mo2_mod_gamma, ask = F)

Checking rhat are equal to one

summary(mo2_mod_gamma)

##  Family: gamma 
##   Links: mu = log; shape = log 
## Formula: MO2 ~ temp_c + order_c + light_dark_c + log_mass_c + salinity_group + (1 | id) 
##          shape ~ salinity_group + order_c
##    Data: slope_tidy_smr (Number of observations: 893) 
##   Draws: 4 chains, each with iter = 8000; warmup = 1000; thin = 2;
##          total post-warmup draws = 14000
## 
## Multilevel Hyperparameters:
## ~id (Number of levels: 58) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.31      0.03     0.26     0.38 1.00     4013     6644
## 
## Regression Coefficients:
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                -2.54      0.06    -2.66    -2.43 1.00     2241
## shape_Intercept           2.58      0.07     2.45     2.72 1.00    12143
## temp_c                   -0.33      0.21    -0.75     0.09 1.00     6431
## order_c                   0.00      0.00    -0.00     0.01 1.00    11445
## light_dark_c              0.13      0.02     0.09     0.18 1.00    12311
## log_mass_c                0.81      0.10     0.61     1.01 1.00     2943
## salinity_group9          -0.10      0.08    -0.26     0.07 1.00     2304
## shape_salinity_group9     0.18      0.10    -0.03     0.37 1.00    12404
## shape_order_c            -0.04      0.01    -0.06    -0.02 1.00    11624
##                       Tail_ESS
## Intercept                 4183
## shape_Intercept          12331
## temp_c                    9409
## order_c                  12155
## light_dark_c             12266
## log_mass_c                5032
## salinity_group9           4394
## shape_salinity_group9    12389
## shape_order_c            12394
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Using leave one out (loo) measure of fit, the model appears to preform well, all, but one Pareto k estimates are good (k < 0.7)

loo(mo2_mod_gamma)

## 
## Computed from 14000 by 893 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   2256.8 38.9
## p_loo        71.6  7.1
## looic     -4513.5 77.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     892   99.9%   765     
##    (0.7, 1]   (bad)        1    0.1%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

Model predictions generally align with the observed data

color_scheme_set("red")
p <- pp_check(mo2_mod_gamma, type = "dens_overlay")
p

10.1.6 📈 Results

We did not see a meaningful difference between the routine metabolic rate for fish from the two salinity treatments.

10.1.6.1 Table S1

Table S1: Fixed effect Estimates (β) and 95% Credible Intervals (95% CI)

model_est <- fixef(mo2_mod_gamma, probs = c(0.025, 0.975)) %>%
    as.data.frame() %>%
    tibble::rownames_to_column(var = "Predictor") %>%
    dplyr::mutate(β = round(Estimate, 3), Q2.5 = round(Q2.5, 3), Q97.5 = round(Q97.5,
        3), `95% CI` = paste0("[", Q2.5, ", ", Q97.5, "]"))

model_est %>%
    dplyr::select(Predictor, "β", "95% CI") %>%
    gt()

Predictor	β	95% CI
Intercept	-2.545	[-2.659, -2.428]
shape_Intercept	2.585	[2.45, 2.715]
temp_c	-0.335	[-0.755, 0.087]
order_c	0.002	[-0.002, 0.005]
light_dark_c	0.134	[0.086, 0.182]
log_mass_c	0.805	[0.609, 1.006]
salinity_group9	-0.096	[-0.259, 0.066]
shape_salinity_group9	0.177	[-0.025, 0.372]
shape_order_c	-0.041	[-0.059, -0.023]

Looking at the marginal mean MO2 for each salinity treatment

em_results <- emmeans(mo2_mod_gamma, ~salinity_group)

em_results_df <- em_results %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), ~round(exp(.), 3)))

em_results_df %>%
    gt()

salinity_group	emmean	lower.HPD	upper.HPD
0	0.079	0.070	0.088
9	0.071	0.064	0.080

These are the estimated mean constrast, we have exponentiated the log estimate to get the fold-change. This means the model predicts that sailianty group 0 has ~12% higher metabolic rate than group 9 (1.12)—but this difference is not credible, because the interval for that difference includes no difference (i.e. a ratio of 1.0).

contrast_results <- contrast(em_results, method = "pairwise")

contrast_results_df <- contrast_results %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), ~round(exp(.), 2)))

contrast_results_df %>%
    gt()

contrast	estimate	lower.HPD	upper.HPD
salinity_group0 - salinity_group9	1.1	0.94	1.29

EEM contrasts for the two light phases, where -0.5 is dark, and 0.5 is light

em_results <- emmeans(mo2_mod_gamma, ~light_dark_c, at = list(light_dark_c = c(0.5,
    -0.5)))

contrast_results <- contrast(em_results, method = "pairwise", .rev = TRUE)

contrast_results_df <- contrast_results %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), ~round(exp(.), 2)))

contrast_results_df %>%
    gt()

contrast	estimate	lower.HPD	upper.HPD
light_dark_c0.5 - (light_dark_c-0.5)	1.14	1.09	1.2

Here we we calculate the model-estimated marginal means for MO2 for a unit increase in mass and transform this onto a fold-change scale. Where log(mass) has be centered, a caule of 0 represents the mean mass of fish, and 1 corresponds to a fish with log(mass) = 1 unit above the mean (i.e. 2.7 times large).

log(mass₂) = log(mass₁)+1

If you exponentiate both sides:

mass₂ = mass₁ x exp(1)

mass₂ ≈ mass₁ x 2.718

em_mass <- emmeans(mo2_mod_gamma, ~log_mass_c, at = list(log_mass_c = c(0, 1)))

# Get the contrast (difference in log MO₂), then back-transform
contrast(em_mass, method = "revpairwise") %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), exp)) %>%
    mutate(across(where(is.numeric), ~round(., 2))) %>%
    gt()

contrast	estimate	lower.HPD	upper.HPD
log_mass_c1 - log_mass_c0	2.23	1.86	2.76

Pulling the emmeans draws for our plot

emmeans_draws_rmr <- mo2_mod_gamma %>%
    emmeans(~salinity_group) %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value), salinity_group = as.character(salinity_group))

emmeans_contrast_draws_rmr <- mo2_mod_gamma %>%
    emmeans(~salinity_group) %>%
    contrast(method = "pairwise") %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value))

10.1.6.2 Figure S7 (Figure 1a)

NOTE: This plot is in the main text of the manuscript as Figure 1a

fig_1a <- ggplot() +
  geom_violin(data = slope_tidy_smr,
              aes(x = salinity_group, y = MO2, fill = salinity_group),
              color = NA, alpha = 0.2) +
  geom_jitter(data = slope_tidy_smr,
              aes(x = salinity_group, y = MO2, fill = salinity_group),
              shape = 21, width = 0.3, size = 1, color = "black", alpha = 0.1) +
  geom_point(data = mean_mo2_salinity,
             aes(x = salinity_group, y = mean_mo2, fill = salinity_group),
             size = 4, alpha = 1, stroke = 2, color = "black", shape = 21,
             position = position_nudge(x = 0.05)) +
  stat_pointinterval(data = emmeans_draws_rmr, 
                     aes(x = salinity_group, y = .value),
                     color = "black", fill = "grey", point_interval = "mean_qi", .width = 0.95, shape = 21,  stroke = 2, point_size = 4, alpha = 1,
                     position = position_nudge(x = -0.05)) +
  scale_y_continuous(limits = c(0.00, 0.35)) +  # 👈 Set y-axis range here
  scale_fill_manual(values = c("#4B5320", "#000080")) +
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme(
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
    panel.background = element_rect(fill = "white", colour = NA),
    plot.background = element_rect(fill = "white", colour = NA)
  ) +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "Routine MO2 (mg O2 g/h)"
  )

fig_1a

Figure 1a: Routine metabolic rate (i.e. MO₂ (mg^-1 O₂ h^-1) measured during SMR meassurments) plotted by salinity treatment. The small transparent points are the observed values, the shaded area behind the points is the a kernel density of the observed data, the large coloured point (to the right) is the observed mean, the large grey point with error bars (to the left) is the model estimated marginal mean (eemean) 95% Credible Intervals (95% CI).

# ggsave(filename = paste0(output_fig_wd, './ms_figures./figure-1a.pdf'), plot
# = fig_1a, width = 210/2, height = 297/2, # Specify the height of the plot
# units = 'mm')

10.2 SMR

10.2.1 Formating and scaling data

Here we are filtering the data frame to have only measure per fish for the SMR estimate

scale_list <- c("temp_mean", "log_mass", "cycles")

smr <- slope_tidy_smr %>%
    dplyr::group_by(id) %>%
    dplyr::reframe(temp_mean = mean(temp), log_mass = log_mass[1], mass = mass[1],
        SMR = SMR[1], salinity_group = salinity_group[1], cycles = length(order)) %>%
    dplyr::mutate(across(all_of(scale_list), ~scale(.x, center = TRUE, scale = FALSE),
        .names = "{.col}_c"))

Summary of predictors used in the model below

smr %>%
    dplyr::reframe(temp_range = paste0(round(min(temp_mean, na.rm = TRUE), 2), "–",
        round(max(temp_mean, na.rm = TRUE), 2)), mass_range = paste0(round(min(mass,
        na.rm = TRUE), 2), "–", round(max(mass, na.rm = TRUE), 2)), cycle_range = paste0(round(min(cycles,
        na.rm = TRUE), 2), "–", round(max(cycles, na.rm = TRUE), 2)), number_fish = length(unique(id))) %>%
    gt()

temp_range	mass_range	cycle_range	number_fish
13.87–14.26	0.21–1.6	10–23	58

10.2.2 Model structure

Here we will use a Bayesian Generalised Linear Mixed Model (GLMM) with a Gamma distribution and a log link Gamma(link = "log"), where the shape parameter (α) is also modelled as a function of the salinity group, shape ~ salinity_group. This models MO2 by salinity during the SMR phase to see if the fish held at different salinities have different SMRs. We have also added a few scaled predictors, that may help describe variation in the data, such as fish mass (g; 0.21–1.6) mean temperature (°C; 13.841–14.277), and the number of cycles over which SMR was estimated (1–23). We allowed the the shape parameter (α) to vary as a function of salinity_group to improve fit.

smr_gamma_bf <- bf(SMR ~ temp_mean_c + cycles_c + log_mass_c + salinity_group, shape ~
    salinity_group, family = Gamma(link = "log"))

10.2.3 Prior selection

These are the default priors for the model. We will use these.

priors_default <- get_prior(smr_gamma_bf, data = smr, family = Gamma(link = "log"))
priors_default %>%
    gt()

prior	class	coef	group	resp	dpar	nlpar	lb	ub	source
	b								default
	b	cycles_c							default
	b	log_mass_c							default
	b	salinity_group9							default
	b	temp_mean_c							default
student_t(3, -2.8, 2.5)	Intercept								default
	b				shape				default
	b	salinity_group9			shape				default
student_t(3, 0, 2.5)	Intercept				shape				default

10.2.4 Run model

Here we run the model, I have hashed this out because I have saved the model for quick reloading.

⏭️ Skip this step if you have downloaded the saved model ‘smr_mod_gamma.rds’

Here we reload the model

smr_mod_gamma <- readRDS(file = paste0(output_mods_wd, "./smr_mod_gamma.rds"))

10.2.5 Model diagnostics

Checking model convergence

color_scheme_set("red")
plot(smr_mod_gamma, ask = F)

Using leave one out (loo) measure of fit, the model appears to preform well, two Pareto k estimates falls outside the good range (0.7, 1]

loo(smr_mod_gamma)

## 
## Computed from 14000 by 58 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo    134.6 11.0
## p_loo        11.0  4.8
## looic      -269.2 22.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 0.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     56    96.6%   4028    
##    (0.7, 1]   (bad)       1     1.7%   <NA>    
##    (1, Inf)   (very bad)  1     1.7%   <NA>    
## See help('pareto-k-diagnostic') for details.

Model predictions generally align with the observed data, but there is a lot of uncertainty around this estimate.

color_scheme_set("red")
p <- pp_check(smr_mod_gamma, type = "dens_overlay")
p

10.2.6 📈 Results

We did not see a meaningful difference between the routine metabolic rate for fish from the two salinity treatments.

10.2.6.1 Table S2

Table S2: Fixed effect Estimates (β) and 95% Credible Intervals (95% CI) from a Bayesian Generalised Linear Mixed Model (GLMM) with a Gamma distribution and a log link

model_est <- fixef(smr_mod_gamma, probs = c(0.025, 0.975)) %>%
    as.data.frame() %>%
    tibble::rownames_to_column(var = "Predictor") %>%
    dplyr::mutate(β = round(Estimate, 3), Q2.5 = round(Q2.5, 3), Q97.5 = round(Q97.5,
        3), `95% CI` = paste0("[", Q2.5, ", ", Q97.5, "]"))

model_est %>%
    dplyr::select(Predictor, "β", "95% CI") %>%
    gt()

Predictor	β	95% CI
Intercept	-2.740	[-2.879, -2.603]
shape_Intercept	2.010	[1.428, 2.532]
temp_mean_c	-0.364	[-1.51, 0.802]
cycles_c	-0.005	[-0.026, 0.018]
log_mass_c	0.816	[0.597, 1.034]
salinity_group9	-0.086	[-0.283, 0.111]
shape_salinity_group9	0.117	[-0.733, 0.945]

Looking at the marginal mean difference between salinity groups

em_results <- emmeans(smr_mod_gamma, ~salinity_group)

em_results_df <- em_results %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), ~round(exp(.), 3)))


em_results_df %>%
    gt()

salinity_group	emmean	lower.HPD	upper.HPD
0	0.065	0.056	0.074
9	0.059	0.052	0.067

contrast_results_df <- contrast_results %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), ~exp(.)))

contrast(em_results, method = "pairwise") %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), ~round(exp(.), 2))) %>%
    gt()

contrast	estimate	lower.HPD	upper.HPD
salinity_group0 - salinity_group9	1.09	0.89	1.32

Here we we calculate the model-estimated marginal means for SMR for a unit increase in mass and transform this onto a fold-change scale. Where log(mass) has be centered, a vaule of 0 represents the mean mass of fish, and 1 corresponds to a fish with log(mass) = 1 unit above the mean (i.e. 2.7 times large).

em_mass <- emmeans(smr_mod_gamma, ~log_mass_c, at = list(log_mass_c = c(0, 1)))

# Get the contrast (difference in log MO₂), then back-transform
contrast(em_mass, method = "revpairwise") %>%
    as.data.frame() %>%
    mutate(across(where(is.numeric), exp)) %>%
    mutate(across(where(is.numeric), ~round(., 2))) %>%
    gt()

contrast	estimate	lower.HPD	upper.HPD
log_mass_c1 - log_mass_c0	2.26	1.84	2.84

Pulling the emmeans draws for our plot

emmeans_draw_smr <- smr_mod_gamma %>%
    emmeans(~salinity_group) %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value), salinity_group = as.character(salinity_group))

emmeans_contrast_draws_smr <- smr_mod_gamma %>%
    emmeans(~salinity_group) %>%
    contrast(method = "pairwise") %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value))

10.2.6.2 Figure S8 (Figure 1b)

This plot is in the main text of the manuscript as Figure 1b

mean_smr_salinity <- smr %>%
  dplyr::group_by(salinity_group) %>%
  dplyr::reframe(mean_SMR = mean(SMR, na.rm = TRUE))


fig_1b <- ggplot() +
  geom_violin(data = smr,
              aes(x = salinity_group, y = SMR, fill = salinity_group),
              color = NA, alpha = 0.2) +
  geom_jitter(data = smr,
              aes(x = salinity_group, y = SMR, fill = salinity_group),
              shape = 21, width = 0.3, size = 1, color = "black", alpha = 0.1) +
  geom_point(data = mean_smr_salinity,
             aes(x = salinity_group, y = mean_SMR, fill = salinity_group),
             size = 4, alpha = 1, stroke = 2, color = "black", shape = 21,
             position = position_nudge(x = 0.05)) +
  stat_pointinterval(data = emmeans_draw_smr,
                     aes(x = salinity_group, y = .value),
                     color = "black", fill = "grey", point_interval = "mean_qi", .width = 0.95, shape = 21,  stroke = 2, point_size = 4, alpha = 1,
                     position = position_nudge(x = -0.05)) +
  scale_y_continuous(limits = c(0.00, 0.35)) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme(
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 1),
    panel.background = element_rect(fill = "white", colour = NA),
    plot.background = element_rect(fill = "white", colour = NA)
  ) +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "Standard metabolic rate (SMR; mg O2 g/h)"
  )

fig_1b

labs( y = expression(“mg O”[2]~“h”^-1), x = “Dissolved oxygen (DO; % saturation)” ) + Figure 1b: The standard metabolic rate estimate (mg^-1 O₂ h^-1) plotted by salinity treatment. The small transparent points are the observed values, the shaded area behind the points is the a kernel density of the observed data, the large coloured point (to the right) is the observed mean, the large grey point with error bars (to the left) is the model estimated marginal mean (eemean) 95% Credible Intervals (95% CI).

# ggsave(filename = paste0(output_fig_wd, './ms_figures./figure-1b.pdf'), plot
# = fig_1b, width = 210/2, height = 297/2, # Specify the height of the plot
# units = 'mm')

10.3 Incremental regression analyses

Here we are following the methods Urbina et al. (2012)^[1] with an incremental regression analyses, in order to determine the best fit for the MO₂ vs O₂ data

This analysis approach evaluates the relative ‘fit’ of each polynomial order equation starting at zero and increasing to the third order, permitting a mathematical assessment of whether fish were oxyconforming or oxyregulatoring. If the data is best fitted/predicted by a single linear relationship (1^st-order polynomial) with a positive slope, this would suggests the fish were oxyconforming. Alternately, if the relationship is best modelled by a flat regression (0^th-order polynomial), or a higher order polynomial (2^nd or 3^rd-order polynomial) the fish is likely oxyregulatoring.

10.3.1 Building Bayesian regressions

Here we are using a Bayesian approach to model fitting with brm. These models take a long time to run, so I have saved them and re-loaded them to save time. I have also saved the summary data produced from the models, to save time, you can simply skip the hashed code and input the resulting summary data.

We will run our custom function, bayes_incremental_regression_by_id.

Let’s make define the output directory

output_mods_bayes_wd <- paste0(wd, "./output-mod./bayes-regs")

⏭️ Skip this step if you have already run this once, or have downloaded the saved models or saved data files from GitHub (that’s why its hashed out). ⌛This code takes a while to run ⌛

During this part of the script I received an error, ‘Caused by error in socketConnection():’. I think the system may be hitting a limit on the number of simultaneous socket connections. If you this part of the code and also get this issue, try reducing the number of parallel workers: plan(multisession, workers = 2), or run again after the error.

# ids <- slope_tidy %>% dplyr::distinct(id) %>% pull(id) plan(multisession)
# future_map( ids, bayes_incremental_regression_by_id, id_name = 'id', data =
# slope_tidy, response = 'MO2_g', predictor = 'DO', seed_number = 143019,
# save_models = TRUE, mod_output_wd = output_mods_bayes_wd ) plan(sequential)

Load all models and store in a list, will use a lot of memory. You can also skip this step and load the resulting data frames below. I am using the custom function load_rds, so we can compare them and generate predictions.

⏭️ Skip this step if you have downloaded the saved data pulled from these models have also hased this out, bayes_reg_mods_fit.csv’ and ‘bayes_reg_mods_predictions.csv’. These are in the mod-data directory.

# bayes_reg_mods <- load_rds(model_dw = output_mods_bayes_wd)

10.3.2 Model fits

⏭️ Skip this step if you have downloaded 💿 bayes_reg_mods_fit.csv.

It gets model fit parameters loo and r2 using the custom function, incremental_regression_bayes_fits. ⌛ This code takes a while to run ⌛

# bayes_reg_mods_fit <- incremental_regression_bayes_fits(models =
# bayes_reg_mods) write.csv(bayes_reg_mods_fit, paste0(mod_data_wd,
# './bayes_reg_mods_fit.csv'), row.names = FALSE)

Reading in this model fit data frame, in the case you did not load in all the models.

bayes_reg_mods_fit <- read.csv(paste0(mod_data_wd, "./bayes_reg_mods_fit.csv"))

Selecting the best fitting model

elpd_loo, or the expected log pointwise predictive density for leave-one-out cross-validation, is a metric used in Bayesian model evaluation to assess the predictive accuracy of a model. The elpd_loo is an approximation of how well the model is expected to predict new data, based on leave-one-out cross-validation. Higher elpd_loo values indicate better predictive performance.

best_fit_bayes_reg <- bayes_reg_mods_fit %>%
    dplyr::group_by(id) %>%
    dplyr::mutate(elpd_loo_rank = rank(-elpd_loo)) %>%
    dplyr::select(id, model_type, elpd_loo, r2, elpd_loo_rank, r2_q2.5, r2_q97.5,
        estimate_DO, conf.low_DO, conf.high_DO) %>%
    ungroup()

10.3.3 Model predictions

⏭️ Skip this step if you have downloaded 💿 bayes_reg_mods_predictions.csv.

It pulls our model predictions using a custom function bayes_mod_predictions.

# bayes_reg_mods_predictions <- bayes_mod_predictions(models = bayes_reg_mods,
# original_data = slope_tidy) write.csv(bayes_reg_mods_predictions,
# paste0(mod_data_wd, './bayes_reg_mods_predictions.csv'), row.names = FALSE)

Reading in the predicted data

bayes_reg_mods_predictions <- read.csv(paste0(mod_data_wd, "./bayes_reg_mods_predictions.csv"))

We are going to combined this with our best fitting model df, so we know how they ranks for LOO.

bayes_reg_mods_predictions <- left_join(bayes_reg_mods_predictions, best_fit_bayes_reg,
    by = c("id", "model_type"))

10.3.4 📈 Results

10.3.4.1 Model selection summary

In most cases, the relationship between MO₂ and DO was best modelled with a 2^nd-order polynomial (n = 22, 38% of fish), followed by a 3^rd-order polynomial (n = 15, 26%), 1^st-order polynomial (n = 11, 19%), and finally 0^th-order polynomial (n = 10, 17%).

For the two most common models, 2^nd- and 3^rd-order polynomials, this could suggest the presence of a critical oxygen threshold (Pcrit) where the relationship between MO₂ and O₂ changes. To confirm the presence of a Pcrit, we need to validated the shape of the polynomials (see Pcrit model below). In any case, these type of models are indicative of oxyregulator.

The next most common model was is 1^st-order polynomial. In the case of the 1^st-order polynomials, it suggest the presences of linear relationship between o2 and MO2, which is indicative of oxyconformer. However, to be true evidence of a oxyconformer this relationship should be positive (i.e. as O₂ falls MO₂ also falls). 13 of the 18 individuals best modelled with a linear function had positive estimates with credible intervals that did not overlap with zero (Table S3).

Lastly, 0^th-order was the least common (n = 3, 5%), it suggests that MO₂ does not show a statistically significant dependence on the O₂. In other words, the metabolic rate does not adjust based on oxygen availability, and there is no clear critical oxygen threshold (Pcrit) where the relationship changes. This is indicative of a oxyregulator.

best_mod <- best_fit_bayes_reg %>%
    dplyr::filter(elpd_loo_rank == 1)

total_fish <- nrow(best_mod)

mod_summary_table <- best_mod %>%
    dplyr::group_by(model_type) %>%
    dplyr::reframe(n = length(id), percent = round((n/total_fish) * 100, 2)) %>%
    dplyr::mutate(best_model_name = case_when(model_type == "lm_0" ~ "0th-order polynomial",
        model_type == "lm_1" ~ "1st-order polynomial", model_type == "lm_2" ~ "2nd-order polynomial",
        model_type == "lm_3" ~ "3rd-order polynomial", TRUE ~ "ERROR"))


table_bwm <- mod_summary_table %>%
    dplyr::select(best_model_name, everything(), -model_type) %>%
    gt() %>%
    cols_align(align = "center", columns = everything())

table_bwm

best_model_name	n	percent
0th-order polynomial	10	17.24
1st-order polynomial	11	18.97
2nd-order polynomial	22	37.93
3rd-order polynomial	15	25.86

Summary of fish best model with a linear function.

10.3.5 Table S3

Table S3: Ten fish that were best modelled with a linear function, showing r², and estimate (). Only two fish had positive estimates with credible intervals that did not overlap with zero, which are highlighted as conforming in the table. Thus in total, 13 of 58 fish showed sufficient evidence to be conforming.

table_lm_1 <- best_mod %>%
    dplyr::filter(model_type == "lm_1") %>%
    dplyr::mutate(r_sq_ci = paste0(format(round(r2, 3), scientific = FALSE), " (",
        format(round(r2_q2.5, 3), scientific = FALSE), " to ", format(round(r2_q97.5,
            3), scientific = FALSE), ")"), est_ci = paste0(format(round(estimate_DO,
        4), scientific = FALSE), " (", format(round(conf.low_DO, 4), scientific = FALSE),
        " to ", format(round(conf.high_DO, 4), scientific = FALSE), ")"), conformer = if_else(conf.low_DO >
        0, "Conforming", "Not conforming")) %>%
    dplyr::select(id, estimate_DO, r_sq_ci, est_ci, conformer) %>%
    dplyr::arrange(conformer, desc(est_ci)) %>%
    dplyr::ungroup()

first_order_ids <- table_lm_1 %>%
    distinct(id) %>%
    pull(.)

mean_mo2 <- slope_tidy %>%
    dplyr::filter(id %in% first_order_ids) %>%
    dplyr::group_by(id) %>%
    dplyr::reframe(mean_MO2_g = mean(MO2_g, na.rm = TRUE))

table_lm_1 <- left_join(table_lm_1, mean_mo2, by = "id")

table_lm_1 <- table_lm_1 %>%
    dplyr::mutate(percent_change = (estimate_DO/mean_MO2_g) * 100) %>%
    dplyr::select(id, r_sq_ci, est_ci, percent_change, conformer)


table_lm_1 %>%
    gt() %>%
    cols_align(align = "center", columns = everything()) %>%
    cols_label(id = "Fish ID", r_sq_ci = "r2 (CI)", est_ci = "Estimate (CI)", percent_change = "Percentage change",
        conformer = "Evidence of oxyconforming")

Fish ID	r2 (CI)	Estimate (CI)	Percentage change	Evidence of oxyconforming
d_9_25nov_3	0.197 (0.012 to 0.394)	0.0011 ( 0.0003 to 0.0019)	0.7437151	Conforming
a_9_22nov_4	0.312 (0.094 to 0.491)	0.0008 ( 0.0004 to 0.0013)	0.5739880	Conforming
b_9_22nov_1	0.231 (0.055 to 0.396)	0.0005 ( 0.0002 to 0.0008)	0.4310066	Conforming
d_9_25nov_2	0.330 (0.102 to 0.507)	-0.0009 (-0.0013 to -0.0004)	-0.7599381	Not conforming
a_0_25nov_1	0.137 (0.002 to 0.331)	-0.0004 (-0.0007 to 0.0000)	-0.2788575	Not conforming
c_0_22nov_3	0.072 (0.000 to 0.256)	0.0005 (-0.0005 to 0.0016)	0.3013443	Not conforming
a_0_24nov_3	0.058 (0.000 to 0.202)	0.0005 (-0.0003 to 0.0012)	0.2832654	Not conforming
c_9_27nov_2	0.068 (0.000 to 0.255)	0.0003 (-0.0003 to 0.0008)	0.1697994	Not conforming
a_9_22nov_1	0.084 (0.000 to 0.266)	0.0003 (-0.0002 to 0.0008)	0.2272965	Not conforming
b_0_25nov_2	0.108 (0.001 to 0.290)	0.0003 ( 0.0000 to 0.0006)	0.2275231	Not conforming
a_0_24nov_1	0.029 (0.000 to 0.134)	0.0002 (-0.0005 to 0.0009)	0.1072286	Not conforming

10.3.5.1 Ploting all models

Data set with all slopes and which model was best

best_fit <- left_join(slope_tidy, best_mod, by = "id")

Saving all regression, and highlighting the model that has the best fit, based on AIC values. This is called combined_reg_plots.pdf

# # Create a list to store the plots plots <- list() model_preds_list <- list()
# incremental_reg_bayes_wd <- paste0(output_fig_wd,
# './incremental_regressions') for (id_i in ids) { # Filter data for the
# current ID df_i <- bayes_reg_mods_predictions %>% dplyr::filter(id == id_i)
# %>% dplyr::mutate(line_size = if_else(elpd_loo_rank == 1, 2, 1), alpha_value
# = if_else(elpd_loo_rank == 1, 1, 0.4)) x_min <- df_i %>% dplyr::reframe(min =
# min(DO), na.rm = TRUE) %>% dplyr::pull(min) y_max <- df_i %>%
# dplyr::reframe(max = max(MO2_g), na.rm = TRUE) %>% dplyr::pull(max)
# best_weighted_model_i <- best_fit_bayes_reg %>% dplyr::filter(id == id_i &
# elpd_loo_rank == 1) poly_i_name <- best_weighted_model_i %>%
# dplyr::mutate(name = case_when( model_type == 'lm_0' ~ '0th-order',
# model_type == 'lm_1' ~ '1st-order', model_type == 'lm_2' ~ '2nd-order',
# model_type == 'lm_3' ~ '3rd-order', TRUE ~ 'ERROR' )) %>% dplyr::pull(name)
# r_i <- best_weighted_model_i %>% dplyr::mutate(r_sq_ci = paste0(round(r2, 3),
# ' (', round(r2_q2.5, 3), '–', round(r2_q97.5, 3), ')')) %>%
# dplyr::pull(r_sq_ci) # Create the plot plot_i <- ggplot() + geom_ribbon(data
# = df_i, aes(x = DO, y = predicted, ymin = pred_lower, ymax = pred_upper, fill
# = model_type), alpha = 0.1) + geom_line(data = df_i, aes(x = DO, y =
# predicted, colour = model_type, linewidth = line_size, alpha = alpha_value))
# + geom_point(data = df_i %>% dplyr::filter(elpd_loo_rank == 1), aes(x = DO, y
# = MO2_g), alpha = 0.6, colour = 'black', size = 2) +
# scale_colour_manual(values = c('red', 'blue', 'green', 'purple'), labels =
# c('0th Order', '1st Order', '2nd Order', '3rd Order')) +
# scale_size_identity() + # Use the size values directly
# scale_alpha_identity(guide = 'none') + # Remove the alpha legend
# annotate('text', x = x_min, y = y_max, label = paste0('Best fit:
# ',poly_i_name, '\n', 'r2 = ', r_i), hjust = 0, vjust = 1, size = 4) + labs(
# title = paste(id_i), x = 'Dissolved oxygen percentage (DO)', y = 'MO2 (o2
# mg/g/h)', colour = 'Model') + theme_classic() + theme(legend.position =
# 'none') # Store the plot plots[[id_i]] <- plot_i print(plot_i) } pdf(file =
# paste0(incremental_reg_bayes_wd, './combined_reg_plots.pdf'), width = 8,
# height = 6) dev.off() # Close the PDF device

Getting a plot for each best fit regression, and overlaying a global model for that model type.

Here we are grouping fish by best fitting model and getting an average trend. I have hashed the code so you don’t need to re-run the models

⏭️ Skip this step if you have downloaded the global models ‘.rds’ in the ‘bayes-regs-global’ folder

# seed_number = 143019 # Set up parallel backend future::plan(multisession,
# workers = 4) # Adjust workers based on system resources # Define model
# formulas and data model_formulas <- list( lm_0 = bf(MO2_g ~ 1, family =
# gaussian()), lm_1 = bf(MO2_g ~ DO, family = gaussian()), lm_2 = bf(MO2_g ~
# poly(DO, 2), family = gaussian()), lm_3 = bf(MO2_g ~ poly(DO, 3), family =
# gaussian()) ) model_data <- list( lm_0 = best_fit %>% filter(model_type ==
# 'lm_0'), lm_1 = best_fit %>% filter(model_type == 'lm_1'), lm_2 = best_fit
# %>% filter(model_type == 'lm_2'), lm_3 = best_fit %>% filter(model_type ==
# 'lm_3') ) # Run models in parallel with future_map() models <-
# future_map(names(model_formulas), ~{ brm( formula = model_formulas[[.x]],
# data = model_data[[.x]], cores = 4, seed = seed_number, save_pars =
# save_pars(all = TRUE), sample_prior = FALSE, silent = TRUE, file =
# paste0(output_mods_bayes_global_wd, './', .x, '_global') ) }) # Assign model
# names to results names(models) <- names(model_formulas) # Stop parallel plan
# after execution future::plan(sequential)

⏭️ Skip this step if you have downloaded the global_models_pred_df.csv

# global_models_preds <- list() for (mode_name in names(bayes_reg_mods_global))
# { mod_i <- bayes_reg_mods_global[[mode_name]] mode_type =
# str_remove(mode_name, '_global') model_predictions_i <-
# as.data.frame(fitted(mod_i, summary = TRUE)) %>% dplyr::mutate(model =
# mode_name, model_type = mode_type) %>% dplyr::rename(pred_lower = Q2.5,
# pred_upper = Q97.5, predicted = Estimate, pred_error = Est.Error) %>%
# dplyr::select(model, everything()) original_data_i <- best_fit %>%
# filter(model_type == mode_type) %>% dplyr::select(-model_type)
# model_predictions_original_i <- cbind(model_predictions_i, original_data_i)
# global_models_preds[[mode_name]] <- model_predictions_original_i }
# global_models_pred_df <- bind_rows(global_models_preds)
# write.csv(global_models_pred_df, paste0(mod_data_wd,
# '/.global_models_pred_df.csv'), row.names = FALSE)

Load data frame

global_models_pred_df <- read.csv(paste0(mod_data_wd, "/.global_models_pred_df.csv"))

10.3.5.2 Figure S9

Best fit regressions with global models based on that polynomial order.

global_models_pred_df <- global_models_pred_df %>%
  dplyr::mutate(model_name = case_when(
      model_type == "lm_0" ~ "0th-order polynomial",
      model_type == "lm_1" ~ "1st-order polynomial",
      model_type == "lm_2" ~ "2nd-order polynomial",
      model_type == "lm_3" ~ "3rd-order polynomial",
      TRUE ~ "ERROR"
    ))

annotation_data <- mod_summary_table %>%
  dplyr::select(model_type, best_model_name, n)

bayes_reg_mods_predictions_best <- bayes_reg_mods_predictions %>% 
  dplyr::filter(elpd_loo_rank == 1) 

fig_2 <- ggplot() +
  geom_line(data = bayes_reg_mods_predictions_best,
            aes(x = DO, y = predicted, color = id), size = 1) +
  # geom_smooth(data = best_fit,
  #           aes(x = DO, y = MO2_g, color = id), size = 1, alpha = 1, se = FALSE) +
  geom_ribbon(data = global_models_pred_df,
              aes(x = DO, ymin = pred_lower
, ymax = pred_upper, group = model),
              fill = "#FC6C85", alpha = 0.2) +  # Shaded confidence intervals
  geom_line(data = global_models_pred_df,
            aes(x = DO, y = predicted), linewidth = 1.5, color = "#FF007F") +
  facet_wrap(~model_type) +
  scale_color_grey(start = 0.1, end = 0.9) +
  labs(
      title = paste("Best fit regressions grouped by polynomial order"),
      x = "Dissolved oxygen percentage (DO)",
      y = "MO2 (O2 mg/g/h)") +
  theme_classic() +
  theme(legend.position = "none") +
  geom_text(data = annotation_data,
            aes(x = -Inf, y = Inf, label = paste0("italic(n) == ", n)),
            hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, parse = TRUE)
fig_2

Figure S9: Individual regression curves showing the relationship between mass-specific oxygen consumption rate (ṀO₂; mg⁻¹ O₂ h⁻¹) and ambient dissolved oxygen (DO; % saturation), grouped by the best-fitting polynomial order (0ᵗʰ to 3ʳᵈ). Each panel represents fish whose ṀO₂–DO relationship was best modelled by that polynomial order, based on leave-one-out cross-validation (LOO-CV). Grey lines represent individual model fits; the bold pink line shows the group-level trend within each polynomial class. Fish were most commonly best fit by a 2ⁿᵈ-order polynomial (n = 22), followed by 3ʳᵈ-order (n = 15), 1ˢᵗ-order (n = 11), and 0ᵗʰ-order (**n = 10). Higher-order fits (2ⁿᵈ and 3ʳᵈ) and 0ᵗʰ-order fits are generally indicative of oxyregulating responses, while 1ˢᵗ-order fits may suggest oxyconforming behaviour when positively sloped.

# ggsave(filename = paste0(output_fig_wd, './ms_figures./figure-2.pdf'), plot =
# fig_2, width = 210, height = 297/2, # Specify the height of the plot units =
# 'mm')

10.4 O₂crit model

We will calculate O₂crit using Chabot method and function calcO2crit.

This function uses the fifth percentile of the MO₂ values observed at dissolved oxygen levels ≥ 80% air saturation as the criterion to assess low MO₂ values. The algorithm then identifies all the MO₂ measurements greater than this minimally acceptable MO₂ value. Within this sub-set, it identifies the MO₂2 measurement made at the lowest DO and thereafter considers this DO as candidate for breakpoint (named pivotDO in the script). A regression is then calculated using observations at DO levels < pivotDO, and a first estimate of O_2crit is calculated as the intersection of this regression line with the horizontal line representing SMR. The script then goes through validation steps to ensure that the slope of the regression is not so low that the line, projected to normoxic DO levels, passes below any MO₂ values observed in normoxia. It also ensures that the intercept is not greater than zero. Corrective measures are taken if such problems are encountered.

lowestMO2 default is the quantile(MO2[DO >= 80], p=0.05). It is used to segment the data and locate the pivotDO.

ids <- slope_tidy %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

pcrit_model_df_list <- list()
pcrit_models <- list()

for (id_i in ids) {

    df_i <- slope_tidy %>%
        dplyr::filter(id == id_i)

    o2crit <- calcO2crit(Data = df_i, SMR = df_i$SMR[1], lowestMO2 = NA, gapLimit = 4,
        max.nb.MO2.for.reg = 7)

    best_mod_type <- best_mod %>%
        dplyr::filter(id == id_i) %>%
        pull(model_type)

    vaule <- o2crit$o2crit
    lowestMO2 = quantile(df_i$MO2[df_i$DO >= 80], p = 0.05)
    SMR <- o2crit$SMR
    nb_mo2_conforming <- o2crit$Nb_MO2_conforming
    r2 <- o2crit$r2
    method <- o2crit$Method
    p <- o2crit$P[1]
    lethal_point <- min(o2crit$lethalPoints)

    pcrit_model_df <- tibble(id = id_i, pcrit_vaule = vaule, pcrit_smr = SMR, pcrit_lowestMO2 = lowestMO2,
        pcrit_nb_mo2_conforming = nb_mo2_conforming, pcrit_lethal_points = lethal_point,
        pcrit_r2 = r2, pcrit_method = method, pcrit_p = p, best_mod_type = best_mod_type)

    pcrit_model_df_list[[id_i]] <- pcrit_model_df

    pcrit_models[[id_i]] <- o2crit

}

pcrit_model_df <- bind_rows(pcrit_model_df_list)

10.4.1 Ploting O2crit

Here’s we are saving all O_2crit models, even those MO₂-O₂ relationships best modelled by 0^th- and 1^st-order polynomials.

# # Create output directory if needed output_fig_pcrit_chabot_wd <-
# file.path(output_fig_wd, 'model_chabot') if
# (!dir.exists(output_fig_pcrit_chabot_wd)) {
# dir.create(output_fig_pcrit_chabot_wd) } ids <- slope_tidy %>%
# dplyr::distinct(id) %>% dplyr::pull() pcrit_chabot_list <- list() # Open a
# single PDF device pdf(file = file.path(output_fig_pcrit_chabot_wd,
# 'combined_chabot_plots.pdf'), width = 8, height = 6) for (id_i in ids) {
# best_mod_type <- best_mod %>% dplyr::filter(id == id_i) %>% pull(model_type)
# r2 <- pcrit_model_df %>% dplyr::filter(id == id_i) %>% dplyr::mutate(pcrit_r2
# = round(pcrit_r2, 3)) %>% dplyr::pull(pcrit_r2) conforming <- pcrit_model_df
# %>% dplyr::filter(id == id_i) %>% dplyr::mutate(pcrit_nb_mo2_conforming =
# round(pcrit_nb_mo2_conforming, 3)) %>% dplyr::pull(pcrit_nb_mo2_conforming) P
# <- pcrit_model_df %>% dplyr::filter(id == id_i) %>% dplyr::mutate(pcrit_p =
# round(pcrit_p, 3)) %>% dplyr::pull(pcrit_p) SMR <- pcrit_model_df %>%
# dplyr::filter(id == id_i) %>% dplyr::mutate(pcrit_smr = round(pcrit_smr, 3))
# %>% dplyr::pull(pcrit_smr) lowestMO2 <- pcrit_model_df %>% dplyr::filter(id
# == id_i) %>% dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
# dplyr::pull(pcrit_lowestMO2) # Generate and render the plot
# plotO2crit(o2critobj = pcrit_models[[id_i]]) # Add a title mtext( text =
# paste0(id_i), side = 3, line = 2, adj = 0, col = 'blue', font = 2, cex = 1.2
# ) mtext( text = paste0('R2 = ', r2, '; p = ', P, '; CP < SMR = ', conforming,
# '; SMR = ', SMR, '; lowestMO2 = ',lowestMO2), side = 3, line = 1, adj = 0,
# col = 'blue', font = 1, cex = 0.8 ) mtext( text = paste0('Best model fit = ',
# best_mod_type), side = 3, line = 0, adj = 0, col = 'red', font = 1, cex = 0.8
# ) } # Close the PDF device *after* the loop dev.off()

Now printing the fish that were fitted with higher order-polynomials in the htlm document

par(mfrow = c(1, 2), mar = c(4, 4, 4, 2))

pcrit_model_meta <- slope_tidy %>%
    dplyr::group_by(id) %>%
    dplyr::slice(1) %>%
    dplyr::ungroup() %>%
    dplyr::left_join(pcrit_model_df, ., by = "id")

higher_order_ids <- pcrit_model_meta %>%
    dplyr::filter(best_mod_type == "lm_2" | best_mod_type == "lm_3") %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

for (id_i in higher_order_ids) {

    comment <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::slice(1) %>%
        dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
        pull(comment)

    r2 <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
        dplyr::pull(pcrit_r2)

    conforming <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
        dplyr::pull(pcrit_nb_mo2_conforming)

    P <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
        dplyr::pull(pcrit_p)

    SMR <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
        dplyr::pull(pcrit_smr)

    lowestMO2 <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
        dplyr::pull(pcrit_lowestMO2)

    # Generate and render the plot
    plotO2crit(o2critobj = pcrit_models[[id_i]])

    # Add a title
    mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
        font = 2, cex = 1.2)

    mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming), side = 3,
        line = 1, adj = 0, col = "blue", font = 0.5, cex = 0.8)

    mtext(text = paste0("SMR = ", SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 0,
        adj = 0, col = "blue", font = 0.5, cex = 0.8)
}

10.4.2 Visual inspection

We conducted visual inspection of all 58 fish to determine if O_2crits were present. We have loaded this as a data frame 💿 o2crit_check

For this visual inspection, we grouped fish into “no”, “maybe” or “yes”, based on the certainty at which we observed a O_2crits. In the case of yes or maybe, we also estimate dissolved oxygen percentage. All authors did this independently.

10.4.2.1 📈 Results

Let’s figure out how consistent the categorisation was, and what the majority response was. There were 35 cases with no disagreements (60.34%), 23 cases with one or more disagreements (39.66%).

disagreement_summary <- o2crit_check %>%
    group_by(id) %>%
    dplyr::reframe(n_disagreement = n_distinct(pcrit_cat) - 1, cat_disagreement = if_else(n_disagreement ==
        0, 0, 1), majority_cat = names(sort(table(pcrit_cat), decreasing = TRUE))[1],
        certainty = max(table(pcrit_cat)/5 * 100))

total_assesment <- o2crit_check %>%
    distinct(id) %>%
    nrow()

disagreement_summary %>%
    group_by(cat_disagreement) %>%
    dplyr::reframe(n = length(id), percent = round(n/total_assesment * 100, 2)) %>%
    gt()

cat_disagreement	n	percent
0	35	60.34
1	23	39.66

10.4.2.2 Figure S10

Here’s a heat map to show decisions for each fish

# Define custom colours for pcrit_cat
category_colours <- c(
  "no" = "#E03C32",
  "maybe" = "#FFD301",
  "yes" = "#7BB662"
)

# Create heatmap
ggplot(o2crit_check, aes(x = expert, y = id, fill = pcrit_cat)) +
  geom_tile(colour = "white") +
  scale_fill_manual(values = category_colours, na.value = "grey90") +
  labs(
    title = "Visual Assessment of pcrit_cat by Expert and ID",
    x = "Expert",
    y = "Fish ID",
    fill = "Category"
  ) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(size = 7), # adjust size as needed
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank()
  )

Expert deviation from the majority assessment ranged from 22.41 to 6.90%.

o2crit_check_sum <- o2crit_check %>%
    dplyr::left_join(., disagreement_summary, by = "id") %>%
    dplyr::mutate(deviated = if_else(pcrit_cat != majority_cat, 1, 0))

o2crit_check_sum %>%
    group_by(expert) %>%
    summarise(deviations = sum(deviated), deviation_rate = round(deviations/58 *
        100, 2)) %>%
    dplyr::arrange(desc(deviations)) %>%
    gt()

expert	deviations	deviation_rate
ech	13	22.41
jmm	8	13.79
llk	6	10.34
tdc	5	8.62
mg	4	6.90

Let’s make a list of those that had majority assigned as yes for an O_2crit. There was 10 total.

o2crit_visual_yes <- disagreement_summary %>%
    dplyr::filter(majority_cat == "yes") %>%
    dplyr::pull(id)

paste0("Based on the majority visual assessment ", length(o2crit_visual_yes), " fish have possible O₂crit.")

## [1] "Based on the majority visual assessment 10 fish have possible O₂crit."

10.4.3 O₂crit numerical rules

We also used a numerical rule to asses if a O_2crit was observed.

We filtered for only cases were at the lowest O₂ value three consecutive MO₂ measures full below our SMR and fifth percentile of the MO2 values observed at dissolved O2 levels > 80%. In the model output these are called nb_mo2_conforming points.

o2crit_numerical_yes <- pcrit_model_df %>%
    dplyr::filter(pcrit_nb_mo2_conforming > 2) %>%
    pull(id)

paste0("Based on this rule there are ", length(o2crit_numerical_yes), " fish with possible O₂crits.")

## [1] "Based on this rule there are 16 fish with possible O₂crits."

10.4.4 O₂crit numbers

We will now check to see which fish have O_2crit based on the numerical and visual inspections.

There are eight fish that both numerically and visually were determined to have a O_2crit. Two fish with a visually confirmed O_2crit, that did not meet the numerical criteria. Eight fish meet the numerical criteria, but did not pass the visual inspection.

shared_elements <- intersect(o2crit_visual_yes, o2crit_numerical_yes)

unique_to_visual <- setdiff(o2crit_visual_yes, o2crit_numerical_yes)

unique_to_numerical <- setdiff(o2crit_numerical_yes, o2crit_visual_yes)

list(Shared = shared_elements, Unique_to_visual = unique_to_visual, Unique_to_numerical = unique_to_numerical)

## $Shared
## [1] "a_0_26nov_4" "b_0_24nov_1" "b_0_25nov_1" "b_0_25nov_3" "b_0_26nov_1"
## [6] "b_9_21nov_1" "b_9_21nov_2" "b_9_21nov_3"
## 
## $Unique_to_visual
## [1] "c_0_21nov_4" "d_0_21nov_2"
## 
## $Unique_to_numerical
## [1] "a_9_21nov_3" "b_0_24nov_2" "b_0_26nov_3" "b_9_22nov_1" "c_0_21nov_2"
## [6] "c_9_24nov_2" "d_0_21nov_3" "d_9_25nov_3"

10.4.4.1 📈 Results

Looking at mean O_2crit estimates.

pop_mass <- pcrit_model_meta %>%
    dplyr::reframe(mean_mass = mean(mass)) %>%
    pull(mean_mass)

pcrit_model_meta %>%
    dplyr::filter(id %in% shared_elements) %>%
    dplyr::mutate(rel_mass = mass/pop_mass) %>%
    dplyr::reframe(n_o2crit = length(pcrit_vaule), mean_o2crit = mean(pcrit_vaule),
        min_o2crit = min(pcrit_vaule), max_o2crit = max(pcrit_vaule), mean_mass = mean(mass),
        mean_rel_mass = mean(rel_mass), sd_mass = sd(mass), sd_rel_mass = sd(rel_mass)) %>%
    gt() %>%
    fmt_number(columns = where(is.numeric), decimals = 2)

n_o2crit	mean_o2crit	min_o2crit	max_o2crit	mean_mass	mean_rel_mass	sd_mass	sd_rel_mass
8.00	24.94	19.10	32.00	0.46	0.85	0.18	0.34

min_dos <- slope_tidy %>%
    dplyr::filter(!(id %in% shared_elements)) %>%
    group_by(id) %>%
    dplyr::reframe(min_do = min(DO))

min_dos %>%
    dplyr::reframe(mean_min_do = round(mean(min_do), 2), sd_min_do = round(sd(min_do),
        2), min_min_do = round(min(min_do), 2)) %>%
    gt()

mean_min_do	sd_min_do	min_min_do
13.12	3.23	8.82

pcrit_model_meta %>%
    dplyr::filter(id %in% shared_elements) %>%
    dplyr::group_by(salinity_group) %>%
    dplyr::reframe(n_o2crit = length(pcrit_vaule), mean_o2crit = mean(pcrit_vaule),
        min_o2crit = min(pcrit_vaule), max_o2crit = max(pcrit_vaule), mean_mass = mean(mass)) %>%
    gt()

salinity_group	n_o2crit	mean_o2crit	min_o2crit	max_o2crit	mean_mass
0	5	25.32	19.1	32.0	0.3540000
9	3	24.30	21.7	28.9	0.6466667

These are the 8 fish that meet numerical criteria and passed visual inspection.

par(mfrow = c(1, 2), mar = c(4, 4, 4, 2))

for (id_i in shared_elements) {

    comment <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::slice(1) %>%
        dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
        pull(comment)

    r2 <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
        dplyr::pull(pcrit_r2)

    conforming <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
        dplyr::pull(pcrit_nb_mo2_conforming)

    P <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
        dplyr::pull(pcrit_p)

    SMR <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
        dplyr::pull(pcrit_smr)

    lowestMO2 <- pcrit_model_meta %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
        dplyr::pull(pcrit_lowestMO2)

    # Generate and render the plot
    plotO2crit(o2critobj = pcrit_models[[id_i]])

    # Add a title
    mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
        font = 2, cex = 1.2)

    mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
        SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
        font = 1, cex = 0.8)
}

Pulling model method and lethal points

ids <- pcrit_model_meta %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

pcrit_mean_plot <- list()

for (id_i in ids) {
    best_mod_type <- best_mod %>%
        dplyr::filter(id == id_i) %>%
        pull(model_type)
    mod <- pcrit_models[[id_i]]
    data <- mod$origData
    min_leathal <- min(mod$lethalPoints)
    method <- mod$Method
    pcrit <- mod$o2crit
    data_slope <- data %>%
        dplyr::mutate(n_row = 1:nrow(.), leathal = if_else(n_row >= min_leathal,
            "nonregulation", "regulation"), method = method, pcrit = pcrit, best_mod_type = best_mod_type)
    pcrit_mean_plot[[id_i]] <- data_slope
}

pcrit_mean_plot_df <- bind_rows(pcrit_mean_plot)

10.4.4.2 Figure S11

Making a O2crit plot for each of the eight fish that past visual and numerical assessment on the same axis.

pcrit_summary <- pcrit_mean_plot_df %>%
    dplyr::filter(id %in% shared_elements) %>%
    group_by(id) %>%
    summarise(SMR = first(SMR), pcrit = first(pcrit), max_DO = max(DO, na.rm = TRUE),
        method = first(method), .groups = "drop")

smooth_segments <- map_dfr(unique(pcrit_summary$id), function(cur_id) {

    sum_info <- pcrit_summary %>%
        filter(id == cur_id)
    SMR_val <- sum_info$SMR
    method_val <- sum_info$method

    df <- pcrit_mean_plot_df %>%
        filter(id == cur_id, leathal == "nonregulation")

    if (nrow(df) < 2)
        return(NULL)

    # Fit the appropriate model
    if (method_val == "through_origin") {
        model <- lm(MO2 ~ DO + 0, data = df)
    } else {
        model <- lm(MO2 ~ DO, data = df)
    }

    coefs <- coef(model)

    if (method_val == "through_origin") {
        slope <- coefs[1]
        DO_start <- 0
        DO_stop <- SMR_val/slope
    } else {
        intercept <- coefs[1]
        slope <- coefs[2]
        # For LS_reg, start at 0 if prediction is already positive, or at the
        # DO where predicted MO2 becomes 0 if not.
        DO_start <- if (intercept < 0)
            -intercept/slope else 0
        DO_stop <- (SMR_val - intercept)/slope
    }

    # IMPORTANT: Do not bound DO_stop by max(df$DO) if you want to extrapolate
    # DO_stop <- min(DO_stop, max(df$DO, na.rm = TRUE)) # <-- Remove or comment
    # out this line

    # Create a sequence from DO_start to DO_stop
    DO_seq <- seq(DO_start, DO_stop, length.out = 100)

    # Calculate the predicted MO2 values for this sequence
    MO2_pred <- if (method_val == "through_origin") {
        slope * DO_seq
    } else {
        intercept + slope * DO_seq
    }

    tibble(id = cur_id, DO = DO_seq, MO2_pred = MO2_pred)
})

Fig_4 <- ggplot() + geom_point(data = pcrit_mean_plot_df %>%
    dplyr::filter(id %in% shared_elements), aes(x = DO, y = MO2, colour = id, alpha = leathal)) +
    scale_alpha_manual(values = c(regulation = 0.2, nonregulation = 0.8)) + guides(alpha = "none") +
    geom_line(data = smooth_segments, aes(x = DO, y = MO2_pred, colour = id), size = 1.2) +
    geom_segment(data = pcrit_summary, aes(x = pcrit, xend = max_DO, y = SMR, yend = SMR)) +
    geom_segment(data = pcrit_summary, aes(x = pcrit, xend = pcrit, y = 0, yend = SMR)) +
    scale_colour_viridis_d(option = "plasma") + theme_classic() + labs(y = expression("mg O"[2] ~
    "h"^-1), x = "Dissolved oxygen (DO; % saturation)") + theme(legend.position = "bottom")

Fig_4

Figure S10: Oxygen consumption rate (ṀO₂; mg⁻¹ O₂ h⁻¹) and ambient dissolved oxygen (DO; % saturation) plotted for each fish that passed both past visual and numerical assessment for O_2crit. The horizontal dashed lines show SMR, the solid sloped lines represent the linear best fit for ṀO values below regulatory failure. The vertical soild lines represent O_2crit (the intersection of the linear best fit for ṀO values below regulatory failure and SMR).

Let’s make some example plots for each major types of O2-MO2 relationships (i.e. 0^th-, 1^st-, and 2^nd/3^rd- order polynomial).

First an example of a 0th-order polynomial fit (c_9_24nov_4).

fig_2a <- bayes_reg_mods_predictions %>%
  dplyr::filter(id == "c_9_24nov_4", model_type == "lm_0") %>% 
  ggplot() +
  geom_point(aes(x = DO, y = MO2), colour = "black", size = 3, alpha = 0.5) +
  geom_segment(
    aes(x = 0, xend = 100, y = 0, yend = SMR[1]), # SMR
    linetype = "dashed", size = 1, colour = "#0F7AAF") +
  geom_segment(
    aes(x = 0, xend = 100, y = SMR[1], yend = SMR[1]), # conform
    linetype = "dashed", size = 1, colour = "#FFA500") +
  theme_few() +
    labs(
    y = expression("mg O"[2]~"h"^-1),
    x = NULL
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

Second an example of a 1st-order polynomial best fit (d_9_25nov_3).

fig_2b <- bayes_reg_mods_predictions %>%
  dplyr::filter(id == "d_9_25nov_3", model_type == "lm_1") %>% 
  ggplot() +
  geom_point(aes(x = DO, y = MO2), colour = "black", size = 3, alpha = 0.5) +
  # geom_line(aes(x = DO, y = predicted * mass), size = 2, colour = "blue") +
  geom_segment(
    aes(x = 0, xend = 100, y = 0, yend = SMR[1]), # SMR
    linetype = "dashed", size = 1, colour = "#0F7AAF") +
  geom_segment(
    aes(x = 0, xend = 100, y = SMR[1], yend = SMR[1]), # conform
    linetype = "dashed", size = 1, colour = "#FFA500") +
  theme_few() +
  labs(
    y = expression("mg O"[2]~"h"^-1),
    x = NULL
  ) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )

Last an example of a fish best fit with a higher order polynomial (2nd-order).

# Define the example ID
pcrit_example <- "b_0_24nov_1"

# 1. Prepare summary data for the selected ID
pcrit_summary_example <- pcrit_mean_plot_df %>%
  filter(id == pcrit_example) %>%
  summarise(
    SMR = first(SMR),
    pcrit = first(pcrit),
    max_DO = max(DO, na.rm = TRUE),
    method = first(method),
    .groups = "drop"
  )

# 2. Subset the original nonregulation data
df <- pcrit_mean_plot_df %>%
  filter(id == pcrit_example, leathal == "nonregulation")

# 3. Fit model and generate predictions
if (nrow(df) >= 2) {
  SMR_val <- pcrit_summary_example$SMR
  method_val <- pcrit_summary_example$method
  
  if (method_val == "through_origin") {
    model <- lm(MO2 ~ DO + 0, data = df)
    slope <- coef(model)[1]
    DO_start <- 0
    DO_stop <- SMR_val / slope
    MO2_pred <- slope * seq(DO_start, DO_stop, length.out = 100)
  } else {
    model <- lm(MO2 ~ DO, data = df)
    intercept <- coef(model)[1]
    slope <- coef(model)[2]
    DO_start <- if (intercept < 0) -intercept / slope else 0
    DO_stop <- (SMR_val - intercept) / slope
    DO_seq <- seq(DO_start, DO_stop, length.out = 100)
    MO2_pred <- intercept + slope * DO_seq
  }
  
  DO_seq <- seq(DO_start, DO_stop, length.out = 100)
  
  smooth_segment <- tibble(
    id = pcrit_example,
    DO = DO_seq,
    MO2_pred = MO2_pred
  )
}

# 4. Plot the figure
fig_2c <- ggplot() +
  geom_point(
    data = pcrit_mean_plot_df %>% filter(id == pcrit_example),
    aes(x = DO, y = MO2, colour = leathal), size = 3, alpha = 0.5
  ) +
  scale_colour_manual(values = c("regulation" = "black", "nonregulation" = "#c30010")) +
  
  # geom_smooth(data = pcrit_mean_plot_df %>% filter(id == pcrit_example),
  #             aes(x = DO, y = MO2), method = "lm", formula = y ~ poly(x, 3), se = FALSE, alpha = 0.8, size = 2, fullrange = TRUE) +
  
  geom_line(
    data = smooth_segment,
    aes(x = DO, y = MO2_pred), size = 2, colour = "#c30010" #failure  slope
  ) +
  
  
  geom_segment(
    data = pcrit_summary_example,
    aes(x = 0, xend = max_DO, y = SMR, yend = SMR), #SMR
    colour = "#FFA500", linetype = "dashed", size = 1) +
  
  geom_segment(
    data = pcrit_summary_example,
    aes(x = 0, xend = max_DO, y = 0, yend = SMR), #confrom
    linetype = "dashed", size = 1, colour = "#0F7AAF") +
  
  geom_segment(
    data = pcrit_summary_example,
    aes(x = pcrit, xend = pcrit, y = 0, yend = SMR), colour = "darkgreen", size = 2 #Pcrit
  ) +
  
  #coord_cartesian(ylim = c(0, 0.16)) +
  theme_few() +
   labs(
    y = expression("mg O"[2]~"h"^-1),
    x = "Dissolved oxygen (DO; % saturation)"
  ) +
  theme(legend.position = "none")

10.4.4.3 Figure S12 (Figure 2)

Here we are showing the out put of the Chabot O_2crit model,

fig_2_bind <- grid.arrange(fig_2a, fig_2b, fig_2c, ncol = 1)  #bind plot

Figure 2: Figure 2. Representative examples of fish that were best modelled by a 0th-order polynomial (A), 1st-order polynomial (B), and 2nd-order polynomial (C). The orange dashed line shows the SMR, while the blue dashed line shows the theoretical oxyconfroming regression. For plot (C), the red sloped line shows the O2 regulation failure regression (based on the rule-based linear regression method of Claireaux and Chabot, 2016), and the solid horizontal green line shows the O2crit estimate (intersection of O2 failure regression and SMR).

# ggsave(filename = paste0(output_fig_wd, './ms_figures./figure-2.pdf'), plot =
# fig_2_bind, width = 210, height = 297, units = 'mm')

11 Comapring to past data

Here, we have recreated a similar plot to that presented in Urbina, Glover, and Forster (2012)^[1] and have extracted the mean level data from Figure 1a using the metaDigitise package in R^[3]. This data is called urbina_et_al_2012. This allows us to compare the differences in Mo2 and the relationship between Mo2 and O₂.

First making a binned data frame to match Urbina, Glover, and Forster (2012) as closely as possible.

min_o2_kpa <- min(slope_tidy$o2_kpa, na.rm = TRUE)
max_o2_kpa <- max(slope_tidy$o2_kpa, na.rm = TRUE)

time_bin_df <- slope_tidy %>%
  mutate(o2_group = cut(o2_kpa, 
                        breaks = seq(min_o2_kpa, max_o2_kpa, length.out = 13), # 11 intervals, so 12 breakpoints
                        labels = paste0("Group ", 1:12), 
                        include.lowest = TRUE)) %>% 
  dplyr::group_by(o2_group) %>% 
  dplyr::reframe(mean_MO2_g = mean(MO2_g)*31.25,
                 mean_o2_kpa = mean(o2_kpa),
                 n = length(MO2_g)*31.25,
                 MO2_g_sd = sd(MO2_g)*31.25,
                 o2_kpa_sd = sd(o2_kpa))

11.1 Figure S13 (Figure 3)

Now the plot with our mean data and the mean data from Urbina, Glover, and Forster (2012)^[1]

n <- slope_tidy %>%
    dplyr::distinct(id) %>%
    nrow(.)

# Existing plot
fig3 <- time_bin_df %>%
    ggplot() + geom_point(data = slope_tidy, aes(y = MO2_g * 31.25, x = o2_kpa),
    size = 2, color = "grey", alpha = 0.3) + geom_point(data = time_bin_df, aes(y = mean_MO2_g,
    x = mean_o2_kpa), size = 3, colour = "#0E4C92", show.legend = FALSE) + geom_errorbar(data = time_bin_df,
    aes(y = mean_MO2_g, x = mean_o2_kpa, ymin = mean_MO2_g - MO2_g_sd, ymax = mean_MO2_g +
        MO2_g_sd), width = 0.15, colour = "#0E4C92") + geom_errorbarh(data = time_bin_df,
    aes(y = mean_MO2_g, x = mean_o2_kpa, xmin = mean_o2_kpa - o2_kpa_sd, xmax = mean_o2_kpa +
        o2_kpa_sd), height = 0.4, colour = "#0E4C92") + geom_point(data = urbina_et_al_2012,
    aes(x = o2_mean, y = mo2_mean), size = 3, shape = 1, fill = "#D21F3C", color = "#D21F3C",
    stroke = 1) + geom_errorbar(data = urbina_et_al_2012, aes(x = o2_mean, ymin = mo2_mean -
    mo2_sd, ymax = mo2_mean + mo2_sd), width = 0.2, colour = "#D21F3C") + geom_errorbarh(data = urbina_et_al_2012,
    aes(y = mo2_mean, xmin = o2_mean - o2_sd, xmax = o2_mean + o2_sd), height = 0.4,
    colour = "#D21F3C") + annotate("text", x = 0, y = 16, label = bquote(atop("Current data (blue), " *
    italic(n) * " = " * .(n), "Urbina data (red), " * italic(n) * " = 67")), hjust = 0,
    vjust = 1, size = 4) + theme(legend.position = "none", panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black",
        fill = NA, linewidth = 1), panel.background = element_rect(fill = "white",
        colour = NA), plot.background = element_rect(fill = "white", colour = NA)) +
    labs(subtitle = NULL, x = expression("PO"[2] ~ "(kPa)"), y = expression("MO"[2] ~
        "(μmol O"[2] ~ "g"^-1 ~ "h"^-1)) + scale_y_continuous(limits = c(0, 16),
    breaks = seq(0, 16, by = 2)) + scale_x_continuous(limits = c(0, 22), breaks = seq(0,
    22, by = 2))

fig3

Figure 3: Mean and standard error of metabolic rate (MO₂ μmol O₂ g^-1 h^-1) and oxygen concentration (PO₂ kPa) using 12 evenly spaced bins over the O₂ range of observed data (blue filled dots). Compared against the mean and standard error reported in Urbina (2012)^[1] (red open dots). The grey dots are the raw observed data form the present study.

# ggsave(filename = paste0(output_fig_wd, './ms_figures./figure-3.pdf'), plot =
# fig3, width = 210, height = 297/2, # Specify the height of the plot units =
# 'mm')

12 📚 References

^[1] Claireaux, G. and Chabot, D. (2016) Responses by fishes to environmental hypoxia: integration through Fry’s concept of aerobic metabolic scope. Journal of Fish Biology https://doi.org/10.1111/jfb.12833

^[2] Urbina MA, Glover CN, and Forster ME, (2012) A novel oxyconforming response in the freshwater fish Galaxias maculatus. Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology. https://doi.org/10.1016/j.cbpa.2011.11.011

^[3] Pick JL, Nakagawa S, and Noble DWA (2018) Reproducible, flexible and high-throughput data extraction from primary literature: The metaDigitise r package. https://doi.org/10.1111/2041-210X.13118

13 💻 Session information

citation()

## 
## To cite R in publications use:
## 
##   R Core Team (2023). R: A language and environment for statistical
##   computing. R Foundation for Statistical Computing, Vienna, Austria.
##   URL https://www.R-project.org/.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {R: A Language and Environment for Statistical Computing},
##     author = {{R Core Team}},
##     organization = {R Foundation for Statistical Computing},
##     address = {Vienna, Austria},
##     year = {2023},
##     url = {https://www.R-project.org/},
##   }
## 
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.

Here is a detailed list of the session information

sessioninfo::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.2.3 (2023-03-15 ucrt)
##  os       Windows 10 x64 (build 19044)
##  system   x86_64, mingw32
##  ui       RTerm
##  language (EN)
##  collate  English_Australia.utf8
##  ctype    English_Australia.utf8
##  tz       Australia/Sydney
##  date     2025-04-02
##  pandoc   3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package           * version  date (UTC) lib source
##  abind               1.4-5    2016-07-21 [1] CRAN (R 4.2.0)
##  arrayhelpers        1.1-0    2020-02-04 [1] CRAN (R 4.2.3)
##  backports           1.4.1    2021-12-13 [1] CRAN (R 4.2.0)
##  bayesplot         * 1.11.1   2024-02-15 [1] CRAN (R 4.2.3)
##  bitops              1.0-7    2021-04-24 [1] CRAN (R 4.2.0)
##  bridgesampling      1.1-2    2021-04-16 [1] CRAN (R 4.2.3)
##  brms              * 2.21.0   2024-03-20 [1] CRAN (R 4.2.3)
##  Brobdingnag         1.2-9    2022-10-19 [1] CRAN (R 4.2.3)
##  broom               1.0.4    2023-03-11 [1] CRAN (R 4.2.3)
##  broom.helpers       1.15.0   2024-04-05 [1] CRAN (R 4.2.3)
##  broom.mixed       * 0.2.9.5  2024-04-01 [1] CRAN (R 4.2.3)
##  bslib               0.7.0    2024-03-29 [1] CRAN (R 4.2.3)
##  cachem              1.0.7    2023-02-24 [1] CRAN (R 4.2.3)
##  caTools             1.18.2   2021-03-28 [1] CRAN (R 4.2.3)
##  cellranger          1.1.0    2016-07-27 [1] CRAN (R 4.2.3)
##  checkmate           2.3.1    2023-12-04 [1] CRAN (R 4.2.3)
##  cli                 3.6.4    2025-02-13 [1] CRAN (R 4.2.3)
##  coda                0.19-4.1 2024-01-31 [1] CRAN (R 4.2.3)
##  codetools           0.2-19   2023-02-01 [2] CRAN (R 4.2.3)
##  colorspace          2.1-0    2023-01-23 [1] CRAN (R 4.2.3)
##  crosstalk           1.2.1    2023-11-23 [1] CRAN (R 4.2.3)
##  curl                5.2.1    2024-03-01 [1] CRAN (R 4.2.3)
##  data.table        * 1.17.0   2025-02-22 [1] CRAN (R 4.2.3)
##  DEoptimR            1.1-3    2023-10-07 [1] CRAN (R 4.2.3)
##  devtools          * 2.4.5    2022-10-11 [1] CRAN (R 4.2.3)
##  digest              0.6.35   2024-03-11 [1] CRAN (R 4.2.3)
##  distributional      0.4.0    2024-02-07 [1] CRAN (R 4.2.3)
##  doParallel          1.0.17   2022-02-07 [1] CRAN (R 4.2.3)
##  dplyr             * 1.1.4    2023-11-17 [1] CRAN (R 4.2.3)
##  ellipsis            0.3.2    2021-04-29 [1] CRAN (R 4.2.3)
##  emmeans           * 1.11.0   2025-03-20 [1] CRAN (R 4.2.3)
##  estimability        1.5.1    2024-05-12 [1] CRAN (R 4.2.3)
##  evaluate            0.24.0   2024-06-10 [1] CRAN (R 4.2.3)
##  fansi               1.0.4    2023-01-22 [1] CRAN (R 4.2.3)
##  farver              2.1.1    2022-07-06 [1] CRAN (R 4.2.3)
##  fastmap             1.1.1    2023-02-24 [1] CRAN (R 4.2.3)
##  forcats           * 1.0.0    2023-01-29 [1] CRAN (R 4.2.3)
##  foreach             1.5.2    2022-02-02 [1] CRAN (R 4.2.3)
##  formatR             1.14     2023-01-17 [1] CRAN (R 4.2.3)
##  fs                  1.6.4    2024-04-25 [1] CRAN (R 4.2.3)
##  furrr               0.3.1    2022-08-15 [1] CRAN (R 4.2.3)
##  future            * 1.33.2   2024-03-26 [1] CRAN (R 4.2.3)
##  generics            0.1.3    2022-07-05 [1] CRAN (R 4.2.3)
##  ggdist              3.3.2    2024-03-05 [1] CRAN (R 4.2.3)
##  ggplot2           * 3.5.1    2024-04-23 [1] CRAN (R 4.2.3)
##  ggthemes          * 5.1.0    2024-02-10 [1] CRAN (R 4.2.3)
##  globals             0.16.3   2024-03-08 [1] CRAN (R 4.2.3)
##  glue                1.8.0    2024-09-30 [1] CRAN (R 4.2.3)
##  gridExtra         * 2.3      2017-09-09 [1] CRAN (R 4.2.3)
##  gsw                 1.2-0    2024-08-19 [1] CRAN (R 4.2.3)
##  gt                * 0.11.0   2024-07-09 [1] CRAN (R 4.2.3)
##  gtable              0.3.5    2024-04-22 [1] CRAN (R 4.2.3)
##  gtsummary         * 1.7.2    2023-07-15 [1] CRAN (R 4.2.3)
##  highr               0.11     2024-05-26 [1] CRAN (R 4.2.3)
##  hms               * 1.1.3    2023-03-21 [1] CRAN (R 4.2.3)
##  htmltools           0.5.8.1  2024-04-04 [1] CRAN (R 4.2.3)
##  htmlwidgets         1.6.4    2023-12-06 [1] CRAN (R 4.2.3)
##  httpuv              1.6.9    2023-02-14 [1] CRAN (R 4.2.3)
##  httr                1.4.7    2023-08-15 [1] CRAN (R 4.2.3)
##  inline              0.3.19   2021-05-31 [1] CRAN (R 4.2.3)
##  insight             1.1.0    2025-03-01 [1] CRAN (R 4.2.3)
##  iterators           1.0.14   2022-02-05 [1] CRAN (R 4.2.3)
##  janitor           * 2.2.0    2023-02-02 [1] CRAN (R 4.2.3)
##  jquerylib           0.1.4    2021-04-26 [1] CRAN (R 4.2.3)
##  jsonlite            1.8.8    2023-12-04 [1] CRAN (R 4.2.3)
##  knitr               1.42     2023-01-25 [1] CRAN (R 4.2.3)
##  labeling            0.4.3    2023-08-29 [1] CRAN (R 4.2.3)
##  later               1.3.0    2021-08-18 [1] CRAN (R 4.2.3)
##  lattice             0.20-45  2021-09-22 [2] CRAN (R 4.2.3)
##  lazyeval            0.2.2    2019-03-15 [1] CRAN (R 4.2.3)
##  lifecycle           1.0.4    2023-11-07 [1] CRAN (R 4.2.3)
##  listenv             0.9.1    2024-01-29 [1] CRAN (R 4.2.3)
##  loo                 2.8.0    2024-07-03 [1] CRAN (R 4.2.3)
##  lubridate         * 1.9.2    2023-02-10 [1] CRAN (R 4.2.3)
##  magrittr            2.0.3    2022-03-30 [1] CRAN (R 4.2.3)
##  marelac             2.1.11   2023-09-25 [1] CRAN (R 4.2.3)
##  marginaleffects   * 0.21.0   2024-06-14 [1] CRAN (R 4.2.3)
##  MASS                7.3-58.2 2023-01-23 [2] CRAN (R 4.2.3)
##  Matrix              1.5-3    2022-11-11 [2] CRAN (R 4.2.3)
##  matrixStats         1.3.0    2024-04-11 [1] CRAN (R 4.2.3)
##  mclust            * 6.1.1    2024-04-29 [1] CRAN (R 4.2.3)
##  measurements        1.5.1    2023-04-28 [1] CRAN (R 4.2.3)
##  memoise             2.0.1    2021-11-26 [1] CRAN (R 4.2.3)
##  mgcv                1.9-1    2023-12-21 [1] CRAN (R 4.2.3)
##  mime                0.12     2021-09-28 [1] CRAN (R 4.2.0)
##  miniUI              0.1.1.1  2018-05-18 [1] CRAN (R 4.2.3)
##  multcomp            1.4-25   2023-06-20 [1] CRAN (R 4.2.3)
##  munsell             0.5.1    2024-04-01 [1] CRAN (R 4.2.3)
##  mvtnorm             1.2-5    2024-05-21 [1] CRAN (R 4.2.3)
##  nlme                3.1-162  2023-01-31 [2] CRAN (R 4.2.3)
##  oce                 1.8-3    2024-08-17 [1] CRAN (R 4.2.3)
##  opdisDownsampling   1.0.1    2024-04-15 [1] CRAN (R 4.2.3)
##  pacman              0.5.1    2019-03-11 [1] CRAN (R 4.2.3)
##  parallelly          1.37.1   2024-02-29 [1] CRAN (R 4.2.3)
##  pbmcapply           1.5.1    2022-04-28 [1] CRAN (R 4.2.0)
##  performance       * 0.12.0   2024-06-08 [1] CRAN (R 4.2.3)
##  pillar              1.9.0    2023-03-22 [1] CRAN (R 4.2.3)
##  pkgbuild            1.4.4    2024-03-17 [1] CRAN (R 4.2.3)
##  pkgconfig           2.0.3    2019-09-22 [1] CRAN (R 4.2.3)
##  pkgload             1.3.4    2024-01-16 [1] CRAN (R 4.2.3)
##  plotly            * 4.10.4   2024-01-13 [1] CRAN (R 4.2.3)
##  plyr                1.8.8    2022-11-11 [1] CRAN (R 4.2.3)
##  posterior           1.6.0    2024-07-03 [1] CRAN (R 4.2.3)
##  pracma              2.4.4    2023-11-10 [1] CRAN (R 4.2.3)
##  profvis             0.3.8    2023-05-02 [1] CRAN (R 4.2.3)
##  promises            1.2.0.1  2021-02-11 [1] CRAN (R 4.2.3)
##  purrr             * 1.0.1    2023-01-10 [1] CRAN (R 4.2.3)
##  qqconf              1.3.2    2023-04-14 [1] CRAN (R 4.2.3)
##  qqplotr           * 0.0.6    2023-01-25 [1] CRAN (R 4.2.3)
##  QuickJSR            1.3.0    2024-07-08 [1] CRAN (R 4.2.3)
##  R6                  2.5.1    2021-08-19 [1] CRAN (R 4.2.3)
##  ragg                1.2.5    2023-01-12 [1] CRAN (R 4.2.3)
##  RColorBrewer        1.1-3    2022-04-03 [1] CRAN (R 4.2.0)
##  Rcpp              * 1.0.10   2023-01-22 [1] CRAN (R 4.2.3)
##  RcppParallel        5.1.8    2024-07-06 [1] CRAN (R 4.2.3)
##  readr             * 2.1.4    2023-02-10 [1] CRAN (R 4.2.3)
##  readxl            * 1.4.3    2023-07-06 [1] CRAN (R 4.2.3)
##  remotes             2.5.0    2024-03-17 [1] CRAN (R 4.2.3)
##  reshape2            1.4.4    2020-04-09 [1] CRAN (R 4.2.3)
##  respirometry      * 2.0.0    2024-07-18 [1] CRAN (R 4.2.3)
##  rlang               1.1.5    2025-01-17 [1] CRAN (R 4.2.3)
##  rmarkdown           2.21     2023-03-26 [1] CRAN (R 4.2.3)
##  robustbase          0.99-4-1 2024-09-27 [1] CRAN (R 4.2.3)
##  rstan             * 2.32.6   2024-03-05 [1] CRAN (R 4.2.3)
##  rstantools          2.4.0    2024-01-31 [1] CRAN (R 4.2.3)
##  rstudioapi          0.14     2022-08-22 [1] CRAN (R 4.2.3)
##  sandwich            3.1-0    2023-12-11 [1] CRAN (R 4.2.3)
##  sass                0.4.9    2024-03-15 [1] CRAN (R 4.2.3)
##  scales              1.3.0    2023-11-28 [1] CRAN (R 4.2.3)
##  seacarb             3.3.3    2024-02-15 [1] CRAN (R 4.2.3)
##  sessioninfo         1.2.2    2021-12-06 [1] CRAN (R 4.2.3)
##  shape               1.4.6.1  2024-02-23 [1] CRAN (R 4.2.3)
##  shiny               1.8.1.1  2024-04-02 [1] CRAN (R 4.2.3)
##  snakecase           0.11.1   2023-08-27 [1] CRAN (R 4.2.3)
##  SolveSAPHE          2.1.0    2021-05-01 [1] CRAN (R 4.2.0)
##  StanHeaders       * 2.32.9   2024-05-29 [1] CRAN (R 4.2.3)
##  stringi             1.7.12   2023-01-11 [1] CRAN (R 4.2.2)
##  stringr           * 1.5.1    2023-11-14 [1] CRAN (R 4.2.3)
##  survival            3.5-3    2023-02-12 [2] CRAN (R 4.2.3)
##  svUnit              1.0.6    2021-04-19 [1] CRAN (R 4.2.3)
##  systemfonts         1.0.4    2022-02-11 [1] CRAN (R 4.2.3)
##  tensorA             0.36.2.1 2023-12-13 [1] CRAN (R 4.2.3)
##  textshaping         0.3.6    2021-10-13 [1] CRAN (R 4.2.3)
##  TH.data             1.1-2    2023-04-17 [1] CRAN (R 4.2.3)
##  tibble            * 3.2.1    2023-03-20 [1] CRAN (R 4.2.3)
##  tidybayes         * 3.0.6    2023-08-12 [1] CRAN (R 4.2.3)
##  tidyr             * 1.3.0    2023-01-24 [1] CRAN (R 4.2.3)
##  tidyselect          1.2.1    2024-03-11 [1] CRAN (R 4.2.3)
##  tidyverse         * 2.0.0    2023-02-22 [1] CRAN (R 4.2.3)
##  timechange          0.2.0    2023-01-11 [1] CRAN (R 4.2.3)
##  twosamples          2.0.1    2023-06-23 [1] CRAN (R 4.2.3)
##  tzdb                0.3.0    2022-03-28 [1] CRAN (R 4.2.3)
##  urlchecker          1.0.1    2021-11-30 [1] CRAN (R 4.2.3)
##  usethis           * 2.2.3    2024-02-19 [1] CRAN (R 4.2.3)
##  utf8                1.2.3    2023-01-31 [1] CRAN (R 4.2.3)
##  V8                  4.4.2    2024-02-15 [1] CRAN (R 4.2.3)
##  vctrs               0.6.5    2023-12-01 [1] CRAN (R 4.2.3)
##  viridisLite         0.4.2    2023-05-02 [1] CRAN (R 4.2.3)
##  withr               3.0.0    2024-01-16 [1] CRAN (R 4.2.3)
##  xfun                0.38     2023-03-24 [1] CRAN (R 4.2.3)
##  xml2                1.3.6    2023-12-04 [1] CRAN (R 4.2.3)
##  xtable              1.8-4    2019-04-21 [1] CRAN (R 4.2.3)
##  yaml                2.3.7    2023-01-23 [1] CRAN (R 4.2.3)
##  zoo                 1.8-12   2023-04-13 [1] CRAN (R 4.2.3)
## 
##  [1] C:/R
##  [2] C:/Program Files/R/R-4.2.3/library
## 
## ──────────────────────────────────────────────────────────────────────────────